Analytical Tools and Method for Modeling Transport Processes in Fluids

ABSTRACT

An improved computer implemented method for modeling transport processes in fluids is disclosed. The method is adapted to model gas flow including dilute gas flow at high Knudsen numbers. Instead of modeling based on using an infinitesimal fluid element of a continuous medium, the method approximates fluid flow as a model gas flow. The method is based on assuming that each of a plurality of particles, which compose the model gas, travels with a probability between any of two points in space occupied by the model gas by following a ballistic trajectory governed by a law of motion in free space and treating each of the plurality of ballistic particles as a property carrier transporting one or more of mass, momentum, and energy between the points of consecutive collisions. The major steps of the method are indicated in the accompanying flow chart. The method also delivers a new basis for prediction of dynamic evolution of the model gas system by considering a pre-established or known dynamic history of the system during a pre-initial period.

TECHNICAL FIELD

The present invention relates to an improvement of the computerimplemented method for modeling transport processes in fluids.

BACKGROUND ART

Computational Fluid Dynamics (CFD) is used mainly for analyticalpredictions of the fluid transport properties, including velocity, heattransfer coefficient, and pressure; still, CFD has many challenges. Onechallenge is that the users must select the appropriate models tocharacterize their specific problem. However, errors due to uncertaintyin the formulation of the model and often intentional simplification ofthe model may diminish the usefulness of the results of thecomputations. Also, errors in the modeling of the fluids are disturbedwith the choice of the governing equations to be solved and models forthe fluid properties. The sources of uncertainty in physical models arelisted by Mehta, U. B. 1996, “Guide to Credible Computer Simulations ofFluid Flows,” AIAA Journal of Propulsion and Power 12, 5, 940-948. (AlsoAIAA Paper 95-2225). Such include 1) the phenomenon is not understood;2) parameters used in the model are known but with some uncertainty; 3)appropriate models are simplified, thus introducing uncertainty; and 4)experimental confirmation of the models is impossible or is incomplete.The general major disadvantage of any existing or available on marketCFD software consists of using an infinitesimal fluid element viewed asa continuous medium, to which fundamental physical principles areapplied. This approach is in a contradiction with the molecular orparticle nature, thus providing a source of significant uncertainty ininterpreting the results of modeling and calculations.

Simulating a system means to solve its mathematical model and to predictits behavior in different situations based on that solution. Continuoussystems can be modeled with ordinary differential equations (ODE) orpartial differential equations (PDE). Solving an ODE or PDE with aninitial value at an initial time, t₀, results in determining thesystem's evolution. Disadvantageously, this approach does not considerthe dynamic behavior of the fluid system in times earlier than initialtime, t₀, which actually affects the behavior of the fluid system intimes following the initial time.

Fluid as a matter comprises molecules. In contemporary fluid dynamics,it is a continuous distribution of matter and is defined as a substancethat deforms continuously under the action of shear stress. Both liquidsand gases are fluids, and all variables describing fluids are consideredcontinuous functions of the spatial coordinates and time. The flow isassumed to be differentiable and continuous, allowing the mass,momentum, and energy balances, which are commonly known as Navier-Stokes(NS) equations, to be expressed as partial differential equations.

Also, the balance equations can only describe fluid flow approximatelymostly because of the inability to determine the functional distributionfor viscous stress forces and the rate of heat addition and viscousdissipation. Thus, although the conservation laws correctly describe thefluid within a control volume located remotely in a fluid flow, they donot take into consideration the molecular nature of the fluid.Specifically, viscosity in gases arises principally from the moleculardiffusion that transports momentum between layers of flow.

As a general limitation, microscopic or bulk analysis allows thetreatment of a collection of free particles as a continuum wherecollisional mean free path, λ_(f), is much lower than the macroscopiclength scale L of the system. However, the classical continuum approachcannot be applied to situations when the length scale of the system isof the same order or less as compared to the collisional mean free path.From the viewpoint of the classical fluid dynamic, the conditionλ_(f)<<L makes good physical sense, which enables the introduction ofthe concept of fluithe d element with volume V having a linear sizelarge in comparison with λ_(f) but small in comparison with L.

In recent years, flow in microsystems has received much attentionbecause of rapid developments in micromachining technology, in whichmost common components are micro fluidic and micro heat transfersystems. To define a transport process, one must understand conditionsat the boundary of the system. At the current state of the art,specifying the boundary conditions is one of the most difficult parts ofthe analysis. In microfluidic systems, the molecular mean free path iscomparable to the system's characteristic length, so that non-continuumor rarefied effects at molecular scale must be considered.

It is well-known from the prior art that “the Navier-Stokes equations,which are the base for the CFD, are only valid as long as therepresentative physical length scale of the system is much larger thanthe mean free path of the molecules that make up the fluid. In thatcase, the fluid is referred to as a continuum. The ratio of the meanfree path, λ_(f), and the representative length scale, L, is called theKnudsen number, Kn=λ_(f)/L. The NS equations are valid for Kn<0.01.”(www.comsol.com/multiphysics/navier-stokes-equations)

The limitation that the NS equations are valid for Kn<0.01 reduces thefields of CFD applications significantly. Specifically, the CFD methodis inapplicable for modeling of dilute or rarefied gases. It is notuseful for simulations microelectronics and nanoelectronics low-pressureprocessing and plasma processing. It does not simulate a gas flowproperly in microchannels in micromechanical systems (MEMS)applications. “In the slip-transition regimes range, that is,0.001<Kn≤10, the CFD scheme, based on the classical Navier-Stokesequations, becomes inappropriate to describe the gas flow behavior inMEMS devices.” (Hindawi Publishing Corporation, Modelling and Simulationin Engineering, Volume 2016, Article ID 7619746, 9 pages,dx.doi.org/101155/2016/7619746)

The other approach involves fluid simulations based on collisions ofparticles positioned at sites of a three-dimensional lattice. It can beperformed effectively only by massively parallel data processors havingcombinational logic for processing collision rules. This method isenormously complicated and restricted because of the limitation of theavailable computational approaches for modeling real physical systems.

Summing noted above, a satisfactory, accurate calculation of momentum ormass or energy transport properties of free particles interacting bycollisions lies in a lack of an appropriate method for describingadequately a collection of free particles traveling in space viaparticle/molecular collisions between gas particles.

In view of the above, a practical universal analytical tool and methodthat allows effective combination/supplement of macroscopic fluiddynamics with microscopic particulate nature of the matter and capablemodel gases in a wide range of gas pressure for high Knudsen numbers isdesired.

SUMMARY OF INVENTION

One general aspect includes analytical tools and a method for modelingfluid flow such as a gas flow in an unstable and stable condition in therange from the continuum flow regime (Kn<0.01) to the free molecularflow regime (Kn>10).

Another general aspect includes describing adequately a collection offree particles traveling in space via particle/molecular collisions forpredicting one or more of distribution patterns of fluid flow propertiessuch as density, mass flow velocity, and temperature.

More specifically, instead of modeling based on using an infinitesimalfluid element that is treated as a continuous medium, which is the basisfor any CFD application, the inventive method aims to approximate fluidflow in a fluid system as a model gas flow in a model gas system, whichis identical to the fluid system.

According to the method, the model gas is designed to facilitate adistant transport of one or more of properties including one or more ofmass, momentum, and energy, which is implemented by way of a pluralityof particles being in a constant state of mostly random motion andinteraction by collisions, where each of the plurality of particles ofthe model gas is assigned to travel with a probability between any oftwo points of a plurality of points in space occupied by the model gasby following a ballistic trajectory governed by a law of motion in freespace. The ballistic trajectory has a starting point in one of aplurality of points of original collisions and an ending point in one ofa plurality of points of ending collisions, and where each of theplurality of particles is adapted to transport a combination of one ormore of properties including one or more of mass, momentum, and energybetween one of the plurality of points of original collisions and one ofthe plurality of points of ending collisions.

The method further including steps of: specifying a geometry model;defining a net rate of property influx per unit volume from the modelgas system in a general point of a plurality of non-moving points inspace occupied by the model gas at a given time, the net rate ofproperty influx per unit volume, which is formed by way of a pluralityof converging ballistic particles, each of the plurality of convergingballistic particles is selected from the plurality of particles by theballistic trajectory having the ending point in the general point of theplurality of non-moving points in space occupied by the model gas at thegiven time, where the given time is greater than or equal to the initialtime; defining a net rate of property efflux per unit volume from thegeneral point at the given time, the net rate of property efflux perunit volume, which is formed by way of a plurality of divergingballistic particles, each of the plurality of diverging ballisticparticles is selected from the plurality of particles by the ballistictrajectory having the starting point in the general point of theplurality of non-moving points at the given time; forming an integralproperty balance equation for each of one or more properties beingtransported by the plurality of particles in each of the plurality ofnon-moving points, where the integral property balance equation isestablished in a way that, in the general point at the given time, aterm of the net rate of property influx per unit volume is equated tothe sum of the first term of temporal rate of a property change per unitvolume and a second term of the net rate of property efflux per unitvolume; and computing the flow of the model gas, where computing theflow of the model gas includes resolving a system of one or more ofintegral property balance equations in each of the plurality ofnon-moving points at the given time. Implementations of the describedtechniques include hardware, a method or process, or computer softwareon a computer-accessible medium.

In the section 2 “Performing flow computation by using the computer forsimulating the fluid flow,” two best modes for carrying out theinvention (Simulation 1 and Simulation 2) are given to explain thefeasibility of the method and to demonstrate significant improvementfrom the prior art.

Simulation 1 is disclosed for incompressible model gas expanding in thespace between two infinite parallel plates in the open channel at theuniform temperature where the top and bottom plates and the model gasare initially at rest. The velocity profile induced in the model gas dueto the sudden motion of the bottom plate is determined by using acomputer program, including computer readable medium having representedin that program codes.

Simulation 2 is disclosed for incompressible model gas expanding in thespace between two infinite parallel plates in the open channel at theuniform temperature in a case of mixed diffuse and specular particlesscatterings from the plates being at rest. The velocity profile inducedin the model gas due to the pressure gradient along the channel isdetermined by using the computer program, including computer readablemedium having represented in that program codes. Note: Diffusescattering is an act of interaction involving property transfer from agas-solid interface to a scattered particle whereas specular scatteringdoes not involve property exchange between the gas-solid interface to ascattered particle.

The most significant improvement consists in demonstrating, first, thecapability to model and compute flows in rarefied gases and, second, inshowing an improvement of time resolution of computations for unstableflows. Specifically, the results of the numerical calculations inSimulation 1 were made for the nitrogen gas flow in a channel confinedby parallel plates spaced at the distance of 1 meter at the uniformtemperature of 300K and pressure of 0.6 Pa (see FIG. 15). FIG. 15 alsodemonstrates the capability of the inventive method to compute thevelocity flow profile at the testing time as short as one millisecond.The results of the numerical calculations in Simulation 2 were made forthe nitrogen gas flow in a channel confined by parallel plates, whichare characterized by the property accommodation coefficient of 0.7 andspaced at the distance of 0.1 meters at the uniform temperature of 300Kand pressure of 1 Pa (see FIG. 16). The simulation demonstratessignificant improvement from the prior art, which in addition to thecapability to model and compute flows in rarefied gases also consists inthe ability to interpret the effect of the property/momentumaccommodation coefficient. Knudsen number in Simulation 1 is aboutKn˜0.05 and, in Simulation 2, is about Kn˜0.1, the values which do notallow using the CFD scheme based on the classical Navier-Stokesequations as it was referenced above.

Moreover, the improved method will be applicable even in values ofKnudsen number equal to or higher than 10, which is a molecular beamregion. Improvement by increasing the operable range of Knudsen numbersis of great importance for the CFD method and will expand horizons ofCFD applications in the industry and in the technological processing.Finally, both Simulation 1 and Simulation 2 demonstrate improvement fromthe prior art, which consists of the capability to produce usefulresults by much lower computational resources.

Contrary, simulations by prior art at such low gas pressure (μ1 Pa) maybe done (if it is possible at all with a time resolution of 1millisecond and the possibility to predict the steps in the velocityprofiles as shown in FIG. 15) only by a massively parallel dataprocessors and are limited by currently underdeveloped computationalapproaches for modeling real physical statistical systems.

The method offers a basis for a new generation of a Computational FluidDynamics analysis processes since modeling of the fluid system evolutioncan be based not only on the established initial conditions at theinitial time for the fluid system but also on the pre-establisheddynamic history of the fluid system at a period preceding the initialtime, thus increasing certainty and accuracy of prediction of the fluidsystem evolution for extended period of time. The method advances a newgeneration of a Computational Fluid Dynamics analysis processes by asignificant decrease of physical modeling errors due to uncertainty inthe formulation of the model and deliberate simplifications of themodel. The method creates potentially much lower discretization errorssince the flow problem can be formulated by direct application ofintegral or integro-differential equations, thus decreasing errorsoccurred from the representation of the governing flow differentialequations and other physical models as algebraic expressions in adiscrete domain of space and time. A new generation of ComputationalFluid Dynamics software products, which are based on the disclosedanalytical tools and method, are anticipated having capability tomodeling gases from the continuum flow regime (Kn<0.01) to the rarefiedand free molecular flow regime (Kn>10), considerably higher accuracy ofprediction, and lower computation cost.

As will further be appreciated by those of skill in the art, theinvention may be embodied as a method, apparatus/system, and computerprogram products. Other systems, methods, and computer program productsaccording to embodiments will be or become apparent to one with skill inthe art upon review of these drawings and detailed description. It isintended that all such additional systems, methods, and computer programproducts be included within this description, be within the scope of thepresent invention, and be protected by the claims.

BRIEF DESCRIPTION OF DRAWINGS

The invention can be best understood by those having ordinary skill inthe art by reference to the following detailed description whenconsidered with the drawings, which are provided as example schematicsand are not meant to limit the invention, where:

FIG. 1 is a general flowchart of the method steps embodying theprinciples of the method.

FIG. 1A is a block diagram of an alternative approach employed to definedifferent components of the net rate of property influx from the modelgas system in a general non-moving point.

FIG. 2 illustrates ballistic traveling of a particle between twoconsecutive collisions in a fluid system composed of the model gasparticles.

FIG. 3 is a block diagram of an alternative approach employed to definemajor steps needed to calculate the net rate of property influx from themodel gas in the general non-moving point at the given time.

FIG. 4 is a schematic view of a one-dimensional configuration forexplaining the ballistic movement of the particle in a field of theexternal force.

FIG. 5 is a schematic view for explaining a method for determining theaverage magnitude of the instant velocity of the traveling particle withrespect to a nearby particle along a ballistic trajectory of thetraveling particle in the one-dimensional configuration.

FIG. 6 is a schematic view of flow geometry and coordinate systemshowing model gas flow confined between two parallel plates.

FIG. 7 is a schematic view of a model gas system confined betweenparallel plates, shown to illustrate a variety of ballistic trajectoriesof particles in the model gas system having gas-solid interfaces withmixed diffuse and specular components of particle scattering.

FIG. 8 is a block diagram of an alternative approach employed to definethe steps needed to calculate the net rate of property influx per unitvolume from a diffuse gas-solid interface in the general non-movingpoint y at the given time t.

FIG. 9 is a block diagram of an alternative approach employed to definemajor steps needed to calculate the net rate of property influx per unitvolume via one or more of mixed diffuse-specular gas-solid interfaces inthe general non-moving point at the given time.

FIG. 10 is a schematic view of the one-dimensional model gas systembetween parallel infinite plates, shown to illustrate a schematic ofdefining diffuse and specular components of the impingement rate on thegas-solid interfaces the plates.

FIG. 11 is a schematic view of the one-dimensional model gas system,shown to illustrate a method for analytical representation of the netrate of property influx, which is formed by converging ballisticparticles originated from diffuse and specular particle scatterings fromconfining parallel plates.

FIG. 12 is a schematic shown to illustrate a method for analyticalrepresentation of the net rate of property efflux from a stationarypoint in one-dimensional space occupied by the model gas.

FIG. 13 is a block diagram of a word equation of a general propertybalance in the model gas system according to the principles of themethod.

FIG. 14 is a flowchart of the flow calculation by the numerical method.

FIG. 15 is a chart showing the solutions providing the velocity profilesacross between the parallel plates at different moments of time and acomparative analytical solution of a steady-state velocity profile.

FIG. 16 is a chart showing the solutions for the stable velocity profileacross between the parallel plates at different numbers of thesequential approximations with mixed diffuse and specular scatteringfrom the plates and a comparative analytical solution.

FIG. 17 is a schematic representation of the special purpose computerfor implementing the method for simulating fluid flow.

DETAILED DESCRIPTION

Reference will now be made to the preferred embodiments of theinvention, examples of which are illustrated in the drawings. Throughoutthe following detailed description, the same reference numerals refer tothe same elements in all figures.

The transport processes involve the exchange of mass, momentum, energy,and other properties exchanged via exchange by particles. In a moreabstract sense, we view all particle interactions and propertyrandomizing as an exchange of property/information. We have alsorecognized that (1) transport phenomena include any situation thatinvolves a net transfer of property/information between particles andrandomizing properties between interacting particles; (2) a randomizedphysical/statistical property by property exchange between interactingparticles can be taken out into the surrounding gas.

It is to be understood that other embodiments may be utilized andstructural changes may be made without departing from the claimedsubject. Specifically, the following specification is limited to thecase of fluid systems having no internal property generation sourcessuch as energy generation sources. However, the method can begeneralized by including in the balance property/energy generation termwith no difficulty, as apprehended by those skilled in the art.

In this specification, all references are made to absolute time, whichis measured equally in all inertial reference frames. Also, whenreferencing to appropriate laws of motion, Newton's First and SecondLaws of motion are applied to the method. However, the method can begeneralized, for example, by applying a different law of motionapplicable to non-Newtonian or relativistic systems, as apprehended bythose skilled in the art.

Also, in these embodiments, the majority of the transport processes aretreated or viewed as transient or unsteady, implying that processvariables change in time, which is the most common state of a system orprocess. However, the method can be reduced to a steady state or in astable state, in which the values of all variables associated with theprocess do not change with time, by eliminating the time dependence invariables and terms associated with the process with no difficulty, asapprehended by those skilled in the art.

In these embodiments, fluid affected by a uniform body force or externalforce field such as gravitational field is analyzed. It is to beunderstood this limitation is made without departing from the claimedsubject.

Some have recognized that the majority of rules and concepts governingfluid flow in three-dimensional and one-dimensional configurations aresimilar. Therefore, for clarity and exemplary, the description dealswith a method of modeling fluid behavior in the one-dimensional fluidsystem, which comes with simulations in a computer. The method is alsoapplicable to two-dimensional and multi-dimensional particle motions.

Embodiments of the present invention will now be described withaccompanying drawings.

1. Analytical Tools and Method for Modeling Fluid Flow

FIG. 1 is a general flowchart of the method steps embodying theprinciples of the method.

The method for modeling transport processes in fluids in an unsteadystate condition includes, for example, the steps of:

approximating fluid flow as model gas flow (S101);

specifying the geometry model (S102);

defining a net rate of property influx per unit volume from the modelgas system in a general point of a plurality of non-moving points inspace occupied by the model gas at a given time, which is formed by wayof a plurality of converging particles, each of the plurality ofconverging particles is selected from the plurality of particles of themodel gas system by a converging ballistic trajectory having the endingpoint in the general point of the plurality of non-moving points inspace at the given time (S103);

defining a net rate of property efflux per unit volume from the generalpoint of the plurality of non-moving points at the given time, which isformed by way of a plurality of diverging ballistic particles, each ofthe plurality is selected from the plurality of particles of the modelgas system by a diverging ballistic trajectory having the starting pointin the general point of the plurality of non-moving points at the giventime (S104);

defining a temporal rate of property change per unit volume in thegeneral point of the plurality of non-moving points at the given time(S105);

forming an integral property balance equation for each of one or moreproperties being transported by the plurality of particles in each ofthe plurality of non-moving points (S106);

computing the flow of the model gas (S107),

where computing the flow of the model gas includes resolving a system ofone or more of integral property balance equations in each of theplurality of non-moving points at the given time.

1.1 the Fluid Model

The fluid model includes the model of the model gas and the model of thegas-solid interface.

1.1.1 the Model of the Model Gas

The step S101 of approximating fluid flow includes approximating fluidflow including a flow of rarefied gases in a fluid system as a model gasflow in a model gas system being identical to the fluid system.

The model gas as any conventional gas is composed of a plurality ofparticles including molecules being in a constant state of mostly randommotion and interaction by collisions. However, the model gas has someproperties that are different from the properties of conventional gases,which in the prior art are characterized by the generalized propertiesof the ideal gas. The major specific properties of the model gasinclude:

each of the plurality of particles of the model gas is assigned totravel with a probability between any of two points of a plurality ofpoints in space occupied by the model gas by following a ballistictrajectory governed by a law of motion in free space, the ballistictrajectory having a starting point in one of a plurality of points oforiginal collisions and an ending point in one of a plurality of pointsof ending collisions, and

each of the plurality of particles is adapted to transport a combinationof one or more of properties including one or more of mass, momentum,and energy between one of the plurality of points of original collisionsand one of the plurality of points of ending collisions.

The model of the model gas further includes:

treating each point of the plurality of points in space occupied by themodel gas as a point of collisions for each of the plurality ofparticles having the ballistic trajectory having the ending point in oneof the plurality of points at the same time;

treating each of a plurality of the points of collisions as either apoint source for each of a plurality of diverging particles or a pointsink for each of a plurality of converging particles; and

treating each of the plurality of particles moving from the point sourceto the point sink as a property carrier created in the point source attime of original collision by obtaining one or more of properties ofspecific values being intrinsic to the model gas surrounding the pointsource, and ended in the point sink at time of ending collision bytransferring one or more of properties of specific values in the pointsink,

where value of property delivered in the point sink by each of theplurality of converging ballistic particles or value of property takenaway from the point source by each of the plurality of divergingballistic particles is evaluated regarding whether value of the propertyis conserved or changed because of aging and whether value of propertycarried by each of the plurality of particles is modified because ofinteraction with an external field, and

where the velocity of each of a plurality of point sources is assignedto be equal to the mass flow velocity of the model gas in each of aplurality of corresponding points of original collisions at time of theoriginal collision.

1.1.2 the Model of the Gas-Solid Interface

The model of the gas-solid interface includes the steps of:

treating each of a plurality of collisions on a gas-solid interface ofthe model gas system, which has resulted in diffuse particle scatteringfrom the gas-solid interface, as an act of interaction involving aproperty transfer from the gas-solid interface to a scattered particle;and

treating each of a plurality of points of diffuse particle scattering onthe gas-solid interface as a heterogeneous point source for each of aplurality of scattered particles,

where the gas-solid interface reveals mixed diffuse and specularparticle scatterings,

where the velocity of each of a plurality of heterogeneous point sourceson the gas-solid interface is assigned to be equal to the velocity ofthe gas-solid interface in each of a plurality of corresponding pointsof diffuse particle scattering at time of diffuse particle scattering,and

where point source strength of each of the plurality of heterogeneouspoint sources on the gas-solid interface is assigned to be directlyproportional to a property accommodation coefficient in each of theplurality of corresponding points of diffuse particle scattering at timeof diffuse particle scattering, the property accommodation coefficientwhich is probability, for an incident particle, to accommodate one ormore of properties intrinsic to the gas-solid interface and to scatterback in the model gas as a diffuse particle, the property accommodationcoefficient being in a range from zero to one. Note: Diffuse scatteringis an act of interaction involving property transfer from a gas-solidinterface to a scattered particle whereas specular scattering does notinvolve property exchange between the gas-solid interface to a scatteredparticle.

1.1.3 Details of the Fluid Model

FIG. 2 illustrates ballistic traveling of the particle between twoconsecutive collisions in a fluid system composed of the model gasparticles. In the example of FIG. 2, a schematic diagram of the modelgas system composed of identical randomly moving particles B andpositioned in the observer's Cartesian coordinate system 200 is shown.For clarity, the schematic presents a three-dimensional configuration.Note that in the specification and accompanying figures, the observer'scoordinate system is designated by index “200.” To locate the positionof the event, the observer needs only to read the space coordinates atthe location of an event. The clocks at any location within a system aresynchronized. The observer, therefore, allocates space coordinates andtime by recording both the space coordinates and time at the clocknearest the event position. In the example of FIG. 2, for clarity withno limitation of the method, the observer's Cartesian coordinate system200 is oriented so y-axis is along the negative direction of an appliedfield of external force 207 which provides, for each of a plurality ofparticles, an acceleration {right arrow over (g)}.

In the specification, claims, and figures herein, the acceleration issymbolized as g in the one-dimensional configuration. Each of theplurality of particles (shown as black disks), particularly a particle201 (shown as transparent disks), travels in space between its twoconsecutive collisions: an original collision, 202, and an endingcollision, 203, by following a trajectory 204 or 205 governed by anappropriate law of motion such as Newton's laws of motion. In a lack ofexternal force, all ballistic trajectories between consecutivecollisions will be just straight lines, respectively, as indicated by atrajectory 206.

Specifically, referring to FIG. 2, particle 201 obtains a set ofproperties of specific values in a point source 202 at time of theoriginal collision, the set of properties including a set of oneproperty of value, Ψ, which includes but not limited to a scalarproperty of value Ψ or a vector property value {right arrow over (Ψ)}(3D), which is intrinsic to the model gas near point source 202 at thetime of the collision, t_(i)′.

In point source 202 moving with velocity being equal to the mass flowvelocity {right arrow over (u)}(t_(i)′, {right arrow over (r)}′) of themodel gas in {right arrow over (r)}′ at time t_(i)′, each of theplurality of particles in the point source obtains, as a propertiescarrier, a specific magnitude of thermal velocity, v_(T)(t_(i)′, {rightarrow over (r)}′), and one or more of model gas properties beingintrinsic to the model gas surrounding the point of the originalcollision at a time of the original collision and may move in anarbitrary direction relatively to the point source. Each of the propertycarriers converging in point sink 203 transfers in the point of theending collision one or more of a set of properties, which in someembodiments includes just one property Ψ_(in), carried by each of theplurality of particles at the time of the ending collision.

1.1.3 the Model of the Geometry Model and the Dynamic History

Next, the step S102 of specifying the geometry model, for example,further includes steps of:

establishing geometry and boundary conditions,

where geometry and boundary conditions are established during a periodfrom a pre-initial time until the given time; and

establishing a dynamic history of the model gas,

where the dynamic history of the model gas, which is in one or more ofinstant distribution patterns of the model gas properties within themodel gas system at a set of different sequential moments in time, isestablished during a pre-initial period from the pre-initial time untilthe initial time, and

where the pre-initial period is set to be not less than a reliableperiod before the initial time.

Some have recognized that to model the model gas flow at the given time,t, which is greater than or equal to the initial time, t₀, a setting ofinitial flow conditions is needed.

This, for example, is done by establishment of the model gas systemgeometry and model gas flow characteristics during some period advancingthe initial time. The initial time, t₀, is associated with the time ofspecific modification of the model gas system in a specific location,which, for example, includes changing of the temperature or/and velocityof a gas-solid interface, and/or application of an external field thatmodifies a property content. It has also been recognized that thespecific modification of the model gas system in the specific locationcan be sensed in a remote point in space at a delayed time, whichdepends on the nature of a transporting carrier. Specifically,modification of the model gas system in the specific location, forexample on a gas-solid interface, will be delivered in the distant pointin space by diffuse ballistic particles and will initiate a modificationof property in the distant point a local initial time according toEquation (1)

t _(i0) =t ₀+φ_(ar)′,  (1)

where φ_(ar)′ is a traveling time along the ballistic trajectory havinga starting point in one of a plurality of points of original collisionsand an ending point in one of a plurality of points of endingcollisions.

The model gas system geometry, for example, is specified by identifyingboth boundaries, which includes the gas-solid interfaces definingboundary conditions and being characterized by specific properties suchas surface temperature distributions and/or dynamics of geometrymodification, and an available set of external fields of certain spatialconfigurations and certain known dynamics of modification, which arebeing a part of the geometry and are prone to affect the propertiescarried by each particle during its free path traveling.

The step of establishing geometry and boundary conditions in someembodiments is utilized:

in the model gas system with no gas-solid interfaces, by providingpredefined dynamics of modification of an available set of externalfields during the period from a pre-initial time until the given time;

in the model gas system bounded by gas-solid interfaces showing eitherpurely diffuse or mixed diffuse and specular particles scattering, byproviding predefined dynamics of modification of both an available setof external fields and boundary conditions on gas-solid interfaces; and

in the model gas system having an open-channel geometry providing modelgas flow passage through a channel having an inlet opening and an outletopening, by providing additional boundary conditions associated with aneffect of surrounding (model gas pressure gradient, temperaturegradient).

The step of establishing a dynamic history of the model gas systemincludes, for example, defining dynamics of the transported propertymodification during particle's ballistic traveling and defining thedynamic history of model gas properties modifications during the modelgas system evolution.

We have recognized that if the property, Ψ₀, acquired by the particle atthe moment of original collision, t_(i)′, or property, Ψ_(0b), acquiredby the particle at the moment of original diffuse scattering from thegas-solid interface, t_(ib)′, may be vulnerable to some internal agingprocess during particle's ballistic traveling, then the value of theproperty upon arriving in point y at the given time t is changed.Internal radioactive decay, in which an unstable atomic nucleus istransformed into a lighter nucleus accompanied by the emission ofparticles or radiation, is an example of the internal aging process. Themass of the particle is decreased during ballistic traveling time.

By the method, for the particle having initial property Ψ₀, the propertyat time t is calculated from Equation (2) given below:

Ψ_(λ)=Ψ₀ e ^(−λφ) ^(i) ,  (2)

where Ψ_(o) is a starting function characterizing property value carriedby each of the plurality of particles in the starting point of theballistic trajectory at a starting time, Ψ_(λ) is an age functioncharacterizing property value in the ending point of the ballistictrajectory, λ is an aging coefficient, which is in the range fromnegative to positive values including zero value, φ_(i) is a gas-gasballistic traveling time defined as follows

φ_(i) =t−t _(i)′.  (3)

A particle incident on the gas-solid interface is scattered backaccording to either a diffuse, or a specular, or a mixeddiffuse-specular model. A gas-solid interface revealing particlescattering according to the mixed diffuse-specular model ischaracterized by a property accommodation coefficient, σ, which istreated as the probability of an incident particle to be “absorbed” intothe surface with further its accommodation on the gas-solid interface byobtaining one or more of properties of specific values being intrinsicto a point of the incidence followed by “evaporation” back and carryingthe obtained properties.

For the particle having initial property, Ψ_(0b), the property, Ψ_(λb),at time t is calculated from Equation (4) given below:

Ψ_(λb)=Ψ_(0b) e ^(−λφ) ^(ib) ,  (4)

where φ_(ib) is solid-gas ballistic traveling time defined as follows

Ψ_(ib) =t−t _(ib)′.  (5)

Notoriously, when λ=0, then the property is conserved during theballistic traveling time. It is to be understood the limitation that λis constant in the equation above is made without departing from theclaimed subject, and those skilled in the art may calculate theproperty, Ψ_(λ) or Ψ_(λb) by using time dependent decay or growthcoefficient with no difficulty.

We have also recognized that the acceleration field g modifies thevelocity of each particle during ballistic traveling. For clarity withno limitation of the invention, acceleration, g, is chosen to be aconstant within the fluid system. An additional value of the propertyΨ_(g)(t_(i)′,t,g), in the velocity modification, v_(g), during theballistic traveling time, φ_(i), for example, is calculated fromEquation (6) given below:

Ψ_(g)(t _(i) ′,t,g)≡v _(g) =g(t−t _(i)′).  (6)

An example of externally coupled force, which acts on gas particlesdirectly through the field during particle's ballistic traveling andcomes with particle displacement toward the force, is a gravitationforce or electrostatic force if the particles are charged.

Under the method, the property delivered in the point sink at time t orthe property, Ψ_(in)(t_(i)′,y′,t,y), carried by each particle at themoment of the convergence in the point sink positioned in generalnon-moving point y at the given time t, includes both an acquiredproperty, which is obtained by each particle in the point of theoriginal collision in the model gas at the time t_(i)′, Ψ₀(t_(i)′,y′),which, sometimes, can be a subject to decay or growth during theballistic traveling time, and the generated property, which isaccumulated or lost by each particle during its ballistic travelingtime, Ψ_(g)(t_(i)′,t,g) is the property, which is accumulated or lostduring the ballistic traveling time because of interaction with anexternal field such as acceleration, g. For clarity, in this embodiment,Ψ_(g)(t_(i)′,t,g) is a well-defined function having arguments associatedwith traveling time φi and with an external well defined field includingacceleration, g.

Symbolically, if the preceding collision is within the model gas space,the discussed above is expressed:

Ψ_(in)(t _(i) ′,y′,t,y)=φ_(λ)(t _(i) ′,y′,λ,φ _(i))+Ψ_(g)(t _(i)′,t,g),  (7)

where Ψ_(g) is a field function characterizing property value of theknown measure, which is changed during traveling time, φ_(i), because ofinteraction with an external field including modification of momentum ina gravitation field Ψ_(λ) is an age function characterizing propertyvalue in an ending point of the ballistic trajectory, which is definedby Equation 2). Analogously, under the method, the propertyΨ_(b_in)(t_(i)′,y′,t,y) carried by each particle at the moment of theconvergence in a point sink positioned in a general non-moving point yat the given time t, includes both the acquired property, which isobtained by each particle in point y_(b) of its preceding diffusescattering from the gas-solid interface at time t_(ib)′,Ψ_(0b)(t_(ib)′,y_(b)), which, sometimes, can be a subject to decay orgrowth during the ballistic traveling time, and the generated property,which is accumulated or lost by each particle during its ballistictraveling time, Ψ_(g)(t_(ib)′,t,g):

Ψ_(b_in)(t _(ib) ′,y _(b) ,t,y)=Ψ_(λb)(t _(ib) ′,y _(b),Δ,φ_(i))+Ψ_(g)(t_(ib) ′,t,g)  (8)

In Equation (8), Ψ_(λb) is defined by Equation (4). In this embodiment,the generated on the gas-solid interface property Ψ_(Ob)(t_(ib)′,y_(b))is a well-defined function if the mentioned gas-solid interface is apart of the external surface bounding the system. The model gas flow ischaracterized microscopically by the group of particles of mass, m,having the number density of particles in the model gas, n, which arerandomly moving and interacting by collisions and having the magnitudeof the random component of velocity or thermal velocity, v_(T), and theparticle's effective collision cross-section, σ_(c). In addition, nearany point of space occupied by the model gas, this group of particles ischaracterized by the vector of mass flow velocity, u.

By the method, each of the model gas properties, Ψ₀(t_(i),y″) in each ofthe plurality of points in the space of the model gas system, is treateddepending on the time of inquiry or the event time either as property,Ψ⁰⁻(t_(i),y″), when event time, t_(i), is advanced to local initialtime, t_(i0), namely, t_(i)<t_(i0), which is specified by establishing aspatial distribution of each of the model gas properties before thelocal initial time t_(i0), or as property, Ψ(t_(i),y″), when the eventtime is greater than or equal to local initial time t_(i0), namely,t_(i)≥t_(i0), which generally has an unknown spatial distribution of theproperty greater than or equal to initial time t_(i0) and is the subjectof determination. Using Iverson brackets, the discussed above isexpressed:

Ψ₀(t _(i) ,y″)=Ψ⁰⁻(t _(i) ,y″)[t _(i) <t _(i0)]+Ψ(t _(i) ,y″)[t _(i) ≥t_(i0)],  (9)

where each of expressions with square brackets is Iverson bracket, whichconverts a statement in brackets into a number that is one if thestatement is satisfied, and is zero otherwise, Ψ⁰⁻ is a pre-establishedproperty value in each of the plurality of points in space of a modelgas system at the starting time ahead of the local initial time, and Ψis a present property value in each of the plurality of points in thespace of a model gas system at the starting time greater than or equalto the local initial time.

1.2 the Model of Property Balance

The model of property balance is based on using analytical tools, whichare constructed by applying the fluid model. Each of the analyticaltools performs a specific and unique function within the model ofproperty balance, which is, in turn, the basis for the computersimulation of the fluid flow.

1.2.1 a Combined Analytical Tool for Computing the Net Rate of PropertyInflux from the Model Gas in the General Non-Moving Point at the GivenTime

In this specification, claims, and figures, a plain symbol is used fornotating a scalar value and a symbol having a rightwards arrow above isused for notating a vector value in the three-dimensional configuration.In addition, by convention in this patent application and to simplifythe notations, the terms relating to the vectors in one-dimensionalconfiguration, which will appear later, are represented by thecorresponding usual letters. For more clarity, each of the vector termsin one-dimensional configuration will be specifically denoted as avector in the tables afterward in the column “Short description” or inthe specification or the claims. Various terms used in the equations arelisted in Table 1.

TABLE 1 Symbol Short description $\frac{\partial\;}{\partial y}$ thedivergence operator in a one-dimensional configuration B_(in) ^(Ψλ)_Fthe net rate of property influx per unit volume from the surroundingmodel gas t the given time t′_(i) the time of the original collision yposition of the ending point of each of a plurality of convergingballistic particles y′ is a position of the starting point of one of theplurality of converging ballistic trajectories v_(T)(t′_(i), y′) theaverage magnitude of thermal velocity in the starting point at the timeof the original collision Z_(V)(t′_(i), y′) rate of collisions per unitvolume in the starting point at the time of the original collisionv(t′_(i), y′, t, y) a velocity vector in the ending point at the giventime, Q_(i)(t, t′_(i)) probability of traveling along the ballistictrajectory Ψ_(in)(t′_(i), y′, t, y a property value delivered in theending point at the given time by one the first plurality of convergingballistic particles yb1 a lower limit of integration over the volume ofspace occupied by the model gas yb2 an upper limit of integration overthe volume of space occupied by the model gas.

FIG. 1A is a block diagram employed to define different components ofthe net rate of property influx from the model gas system in a generalnon-moving point.

The step S103 of defining the net rate of property influx per unitvolume from the model gas system further includes a step of defining anet rate of total property influx per unit volume from a surroundingmodel gas (S103A),

where the net rate of total property influx per unit volume is formed byway of a first plurality of converging ballistic particles, each of thefirst plurality of converging ballistic particles is selected from theplurality of converging ballistic particles by the ballistic trajectoryhaving the starting point in one of a plurality of points of originalcollisions within space occupied by the model gas, and

where each of the plurality of points of original collisions is treatedas a point source.

By the method, each particle traveling from the point of an originalcollision within the space occupied by the model gas to a destinationpoint at the given time, t, is treated as property carrier delivering inthis point some property obtained in the location and time of theoriginal collision. Again, during the ballistic traveling time, thecarried by the particle property may, however, be vulnerable to someinternal aging process and be changed because of interaction with theexternal field.

It is recognized and is explained afterward that there exists acombination of a specific direction of an initial instant vector ofthermal velocity v_(T)(t_(i)′,y′) and a vector of mass flow velocityu(t_(i)′,y′), which provides an opportunity for each of the selectedparticles to arrive in a given non-moving point at a given time. Theseparticles, which are from original collisions within the entire modelgas system, form a converging flux in the given non-moving point of thespace occupied by the model gas at the given time.

In these embodiments, for certainty, the magnitude of thermal velocityis chosen to be higher than the magnitude of mass flow velocity, namely,v_(T)>|u|. This limitation is made without departing from the claimedsubject.

FIG. 3 is a block diagram of an alternative approach employed to definemajor steps of step S103A needed to calculate the net rate of propertyinflux from the model gas in the general non-moving point y at the giventime t.

The step S103A of defining the net rate of total property influx perunit volume from the surrounding model gas includes steps of:

identifying, for each of the first plurality of converging ballisticparticles, the ballistic trajectory (S103A1);

defining the probability of traveling along the ballistic trajectory(S103A2);

defining a net rate of particle efflux per unit volume from one of theplurality of points of original collisions at the time of the originalcollision (S103A3),

where the ballistic trajectory identifies the time of the originalcollision for each of the plurality of converging ballistic particles,and

where the velocity of the point source for each of the first pluralityof converging ballistic particles is assigned to be equal to the massflow velocity of the model gas in a corresponding point of originalcollisions at time of the original collision;

defining a net rate of property flux per unit area in the general pointat the given time from one of a plurality of point sources (S103A4);

defining a net rate of total property flux per unit area in the generalpoint at the given time from the plurality of point sources (S103A5);and

defining the net rate of total property influx per unit volume in thegeneral point at the given time by applying the divergence operator tothe net rate of total property flux (S103A6).

The step of identifying the converging ballistic trajectory includes:

determining one or more of characteristics of the converging ballistictrajectory including the time of the original collision, an appropriateinstant unit vector defining the ballistic trajectory in starting point,and a velocity vector in the ending point of the ballistic trajectory atthe given time; and

defining a size of an expansion zone.

FIG. 4 is an example of a schematic view of a one-dimensionalconfiguration for explaining a ballistic movement of the particle in afield of the external force. In the example of FIG. 4, the observer'sCartesian coordinate system 200 is oriented with y-axis directed alongthe negative direction of the applied acceleration field, 400. FIG. 4shows a particle 406 having a collision in the event point y″≡y′ at theevent time t_(i) ≡t_(i)′, which resulted in acquiring by the particlethermal velocity, v_(T)(t_(i)′,y′) and mass flow velocity vectoru(t_(i)′,y′), which are intrinsic to the nearby volume of the model gassurrounding point y′ at time t_(i)′.

1.2.1.1 Analytical Tool for Computing the Trajectory and Velocity of aBallistic Movement of the Particle Converging in the General Non-Movingat the Given Time

In one embodiment, the step of identifying a trajectory of a ballisticmovement of the particle that converges in the general non-moving pointy at the given time t includes a step of formulating position vector{tilde over (y)} of particle 406 on trajectory 404 or particle 407 ontrajectory 405 and velocity vector {tilde over (v)} (not shown) at time{tilde over (t)} by applying Equations (10) and (11), respectively,given below:

{tilde over (y)}=y({tilde over (t)})=y′+[v _(T)(t _(i) ′,y′)n _(i) +u(t_(i) ′,y′)]({tilde over (t)}−t _(i)′)+½g({tilde over (t)}−t_(i)′)²  (10)

and

{tilde over (v)}=v({tilde over (t)})=v _(T)(t _(i) ′,y′)n _(i) +u(t _(i)′,y′)+g({tilde over (t)}−t _(i)′).  (11)

where t≥{tilde over (t)}≥t_(i)′. In the equations above, n_(i)=±1 is aunit vector of arbitrary direction.

The step of determining the time of the original collision, t_(i)′, inpoint y′ for each of the ballistic traveling particles includes a stepof resolving the equation of motion, Equation (12) given below, withrespect to the ballistic traveling time, φ_(i):

½gφ _(i) ² +v _(T)(t _(i) ′,y′)n _(i) +u(t _(i) ′,y′)φ_(i)+y′−y=0,  (12)

which is obtained by assigning {tilde over (t)}=t and substitutionφ_(i)=t−t_(i)′ in Equation (10) and rearrangement of the terms.

The step of defining an appropriate instant unit vector directingthermal velocity component of each particle in point y′ at time t_(i)′includes steps of:

presenting Equation (12) in the following form:

φ_(i) v _(T) n _(i) =y−y _(i) ^(c),  (13)

where

y _(i) ^(c) =y′+φ _(i) u(t _(i) ′,y′)+½g(φ_(i))²;  (14)

deriving n_(i) from Equation (13), which is given as

$\begin{matrix}{{n_{i} = \frac{y - y_{i}^{c}}{{y - y_{i}^{c}}}},} & (15)\end{matrix}$

where referring to FIG. 4, n_(i) is interpreted as a unit exteriorfacing outside of control volume or an expansion zone 401 having a pointof origin y_(i) ^(c) at time t. In this, vector y_(i) ^(c) isinterpreted as a position vector of the center of the expansion zone attime t or the position of a particle having zero magnitude of thermalvelocity in point y′ at time t_(i)′, which is observed at time t.

The step of defining the velocity of each particle reaching the generalnon-moving point y at the given time t includes a step of applyingEquation (16) given below:

v(t,y)=v _(T)(t _(i) ′,y′)n _(i) +u(t _(i) ′,y′)+gφ _(i),  (16)

which is obtained by assigning t=t and substitution of

φ_(i) =t−t _(i)′  (17)

in Equation (11) given above and rearrangement of the terms.

The step of defining the size of the expansion zone includes a step ofexecuting multiplication of Equation (13) on itself resulted andaveraging as:

φ_(i) ² v _(T) ²=(y−y _(i) ^(c)),  (18)

where the size or half width of the expansion zone 401 at time t iscomputed as follows

Y _(i) ^(exp) =v _(T)(t _(i) ′,y′)φ_(i) =|y−y _(i) ^(c)|  (19)

and the velocity of the center of the expansion zone, v_(i) ^(c), 401 iscomputed as

v _(i) ^(c) =u(t _(i) ′,y′)+gφ _(i).  (20)

1.2.1.2 Analytical Tool for Computing the Probability of BallisticTraveling Along the Ballistic Trajectory of the Particle Converging inthe General Non-Moving at the Given Time

The step of formulating the probability of ballistic traveling along theballistic trajectory includes steps of:

defining the average magnitude of the relative velocity or the velocityof the traveling particle with respect to one of nearby particles at aparticular point of a trajectory; and

representing the probability of traveling along the ballistic trajectoryfrom the starting point to the ending point in the one-dimensionalconfiguration by following Equation (21) given below:

Q _(i)(t,t _(i)′)=exp(−∫_(t) _(i) _(′) ^(t) P _(c)({tilde over(y)}({tilde over (t)}))v _(rel)({tilde over (y)}({tilde over(t)}))d{tilde over (t)}),  (21)

where {tilde over (t)} is a parametric time t_(i)′<{tilde over (t)}≤t,{tilde over (y)}({tilde over (t)}) is a trajectory point of theballistic trajectory at the parametric time {tilde over (t)},v_(rel)({tilde over (y)}({tilde over (t)})) is an average magnitude ofrelative velocity in the trajectory point, and P_(c)({tilde over(y)}({tilde over (t)})) is average number collisions per unit length inthe trajectory point, and

where Equation (21) is obtained by

defining Q({tilde over (t)}+d{tilde over (t)}) and Q({tilde over (t)})as probabilities that a given particle traveling along the ballistictrajectory has survived at times {tilde over (t)}+d{tilde over (t)} and{tilde over (t)}, respectively;

defining y({tilde over (t)}) as the position of the converging ballisticparticle along the ballistic trajectory designated at time {tilde over(t)}, and t_(i)′ is the time of the original collision,

representing P_(c)(y({tilde over (t)},t_(i)′)) by Equation (22) givenbelow:

P _(c)({tilde over (y)}({tilde over (t)}))=σ_(c) n({tilde over(y)}({tilde over (t)})),  (22)

where σ_(c) is the cross-section of collisions and n({tilde over(y)}({tilde over (t)})) is particle density in point {tilde over(y)}({tilde over (t)}) of the converging ballistic trajectory at time{tilde over (t)};

defining v_(rel) as an average magnitude of the velocity of thetraveling particle with respect to a nearby particle along the ballistictrajectory of the traveling particle in a particular point of thetrajectory at the designated time {tilde over (t)};

establishing the following quantity:

Q({tilde over (t)}+d{tilde over (t)})=Q({tilde over (t)})[1−P_(c)(y({tilde over (t)}))v _(rel)(y({tilde over (t)}))d{tilde over(t)}];  (23)

solving a differential based Equation (23) given above and applyinglimits for integration from t_(i)′ to t.

1.2.1.3 Analytical Tool for Computing the Average Magnitude of theRelative Velocity of the Particle Traveling Along the BallisticTrajectory

FIG. 5 is a schematic view for explaining a method for determining theaverage magnitude of the instant velocity of the traveling particle withrespect to a nearby particle along the ballistic trajectory of thetraveling particle in the one-dimensional configuration.

The embodiments of this and the following sub-sections dealing with theaverage magnitude of the relative velocity or the velocity of thetraveling particle with respect to on of nearby particles are explainedwhere traveling particles are non-relativistic. The present invention,however, is not limited to the case but can be generalized to thetransport processes involving relativistic particles, as would beapparent to one skilled in the art.

The step of defining the average magnitude of relative velocity oraverage velocity of a traveling particle 505 with respect to one ofnearby particles 506 at a particular point of the trajectory {tilde over(y)} at specified time {tilde over (t)} is calculated by steps of:

defining an instant magnitude of the velocity of one of the plurality ofconverging ballistic particles in the trajectory point with respect toone of nearby passed particles in the trajectory point, which includes:

defining in a rest frame the instant velocity vector of the travelingparticle 505 in a particular point of the trajectory y({tilde over (t)})at specified time {tilde over (t)}, v₁, by Equation (24) given below:

v ₁ =v _(T)(t _(i) ′,y′)n ₁ +u(t _(i) ′,y′)+g({tilde over (t)}−t_(i)′),  (24)

where n₁ is a unit vector of a specific direction having the point oforigin y′;

specifying instant vector-velocity in the rest frame, v₂({tilde over(t)},{tilde over (y)}), of the nearby particle 506 being passed by thetraveling particle 505 at the same time {tilde over (t)} by Equation(25) given below:

v ₂ =n _(i) v _(T)({tilde over (y)}({tilde over (t)}))+u({tilde over(y)}({tilde over (t)})),  (25)

where v_(T)({tilde over (y)}({tilde over (t)})) and u({tilde over(y)}({tilde over (t)})) are thermal and mass flow velocities in the restframe of the model gas in point {tilde over (y)} at time {tilde over(t)}, which are acquired by the nearby particle 506, and n_(i) isdefined as a unit vector of an arbitrary direction having the point oforigin in point {tilde over (y)} at time {tilde over (t)};

forming an instant vector of relative velocity, v_(rl), which isrepresented in the example of FIG. 5 as v_(rl1) or v_(rl2) by connectingthe end of the instant vector-velocity, v₂({tilde over (t)},{tilde over(y)}), positioned in point 501 or 502 and the end of the vector-velocityv₁(t_(i)′,y′,{tilde over (t)}) positioned in point 504, which isformulated by Equation (26) given below:

v _(rl)({tilde over (t)},{tilde over (y)})=v ₁(t _(i) ′,y′,{tilde over(t)})−v ₂({tilde over (y)}({tilde over (t)})),  (26)

where

v _(r) =v ₁(t _(i) ′,y′,{tilde over (t)})−u({tilde over (y)}({tilde over(t)}))  (27)

is the vector-velocity v_(r) shown in FIG. 5 as 503;

defining the magnitude of the instant relative velocity between thetraveling particle 505 and the nearby particle 506 by executing thesquare root of the scalar product of the velocity vector of Equation(26) with itself, which is expressed in formula form by Equation (28)given below:

|v _(rl)|=√{square root over (v _(rl) v _(rl))}=√{square root over (v_(r) ²+[v _(T)({tilde over (y)}({tilde over (t)}))]²−2v _(r) v_(r)({tilde over (y)}({tilde over (t)}))n _(i))},  (28)

where n_(i) is a unit vector of an arbitrary direction, namely, itaccepts value +1 or −1, and where, in the example of FIG. 5, |v_(rl)| isrepresented by v_(rl1) and v_(rl2); and

calculating the average magnitude of relative velocity of the particletraveling along the ballistic trajectory over all possible directions ofthe random instant relative velocity vector component of a passed nearbyparticle by applying Equation (29) given below:

$\begin{matrix}{{v_{rel} = \frac{{{v_{r} - v_{T}}} + {{v_{r} + v_{T}}}}{2}},} & (29)\end{matrix}$

which is obtained by averaging the instant magnitudes of two componentsof a relative velocity between the particles as of Equation (28) sincen_(i) accepts either value +1 or −1.

In most realistic situations, the magnitude of a mass flow velocitycomponent of either traveling particle 505 or the magnitude of a massflow velocity component of nearby particle 506 particle is insignificantin comparison with the magnitude of the thermal velocity of eitherparticle 505 or particle 506.

In a one-dimensional configuration, for example, an approximated averagemagnitude of the relative velocity of the traveling particle at aparticular point of a trajectory, {tilde over (y)}, at the time, {tildeover (t)}, is calculated from Equation (30) given below, which isobtained by substitution of |v_(r)|≅v_(T)(t_(i)′,y′) in Equation (29):

$\begin{matrix}{{{v_{rel}\left( {\overset{\sim}{y}\left( \overset{\sim}{t} \right)} \right)} = \frac{{{{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)} - {v_{T}\left( {\overset{\sim}{y}\left( \overset{\sim}{t} \right)} \right)}}} + \left( {{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)} + {v_{T}\left( {\overset{\sim}{y}\left( \overset{\sim}{t} \right)} \right)}} \right)}{2}},} & (30)\end{matrix}$

where v_(rel)({tilde over (y)}({tilde over (t)})) is the averagemagnitude of relative velocity in the trajectory point, v_(T)(t_(i)′,y′)is the average magnitude of the thermal velocity component in thestarting point of the ballistic trajectory, and v_(T)({tilde over(y)}({tilde over (t)})) is the average magnitude of the thermal velocitycomponent of one of nearby particles in the trajectory point.

We have also recognized that, usually, the magnitudes of the thermalvelocity of nearby particles are approximately identical. Under themethod, in a one-dimensional configuration, the approximated averagemagnitude of relative velocity between each particle moving in arbitrarydirections, which is originated from the proximity to each other, iscalculated from Equation (31) given below, which is obtained bysubstitution of v_(T)({tilde over (t)},{tilde over(y)})=v_(T)(t_(i)′,y′) in Equation (30):

v _(rel)(t _(i) ′,y′)=v _(T)(t _(i) ′,y′).  (31)

1.2.1.4 Analytical Tool for Computing the Net Rate of Property InfluxPer Unit Volume in the General Non-Moving at the Given Time

A method for determining particle and property fluxes by ballisticparticles retreated from a remotely located auxiliary control volume ina field of external forces is further illustrated in FIG. 4. Asexplained before, the model gas is characterized by the particledensity, n, each particle diverges from a point of the originalcollision and converges in a point of the consecutive ending collision.Under the method, upon recognizing that each traveling particle is beingin both diverging and converging conditions simultaneously, the particledensity of diverging ballistic particles in the vicinity of point y′ attime t_(i)′ is chosen to be n/2.

The step of defining the particle flux production rate, Z_(V), or thenet rate of particle efflux per unit volume from a point sourcepositioned in a point of the original collisions and moving with themass flow velocity of the model gas in that point includes, for example,a step of representing Z_(V) by applying Equation (32) given below:

Z _(V)(t _(i) ′,y′)=½P _(c)(t _(i) ′,y′)v _(rel)(t _(i) ′,y′)n(t _(i)′,y′),  (32)

which is obtained by steps of:

defining the vector field of the particle flux, J_(y′) ^(N), along eachof the ballistic trajectories in a point of the space around the pointof the original collision includes a step of representing J_(y′) ^(N) byapplying Equation (33) given below:

J _(y′) ^(N)=½n(t _(i) ′,y′)Q _(i)(t,t _(i)′)v(t _(i) ′,y′,t,y),  (33)

where v is defined by Equation (16), Q_(i) is a survival probabilitydefined by Equation (21), and n(t_(i)′,y′) is particle density at aspecific point y′ at time t_(i)′;

representing, in a coordinate system associated with point y_(i) ^(c)moving with velocity v_(i) ^(c), the vector field of the particle flux,J_(CV) ^(N), through the in a control surface 402 or 403 of FIG. 4 inthe following form:

J _(CS) ^(N)=½n(t _(i) ′,y′)Q _(i)(t,t _(i)′)[v(t _(i) ′,y′)−v _(i)^(c)]=½n(t _(i) ′,y′)Q _(i)(t,t _(i)′)v _(T)(t _(i) ′,y′)n _(i);  (34)

and applying and executing

$\frac{\partial}{\partial y^{\prime}}$

operator to the one-dimensional vector field of Equation (34) followedby shrinking the volume of the auxiliary control volume to an infinitelysmall volume, which, in formula form, is expressed as:

$\begin{matrix}{{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)} = {\left\{ {\frac{\partial}{\partial y^{\prime}}\left\lbrack J_{CS}^{N} \right\rbrack} \right\}_{y\rightarrow y^{\prime}} = {\frac{1}{2}{\left\{ {\frac{\partial}{\partial y^{\prime}}\left\lbrack {{n\left( {t_{i}^{\prime},y^{\prime}} \right)}{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}n_{i}} \right\rbrack} \right\}_{y\rightarrow y^{\prime}}.}}}} & (35)\end{matrix}$

The above is obtained by recognizing that auxiliary control volume 401is confined by inflated control surface 402 and control surface 403 isisolated (see FIG. 4) therefore the total particle efflux throughsurface 402 and surface 403 at time t should be equal to the totalparticle efflux from the point at y′ at time t_(i)′.

The step of formulating the net rate of the property flux originatedfrom original collisions within entire space occupied by the model gasin a point sink positioned in a general non-moving point of the spaceoccupied by the model gas at the given time includes steps of:

representing the property vector flux, J_(y′→y) ^(Ψ)(t,y), originatedfrom original collisions in point y′ at time t_(i)′ and being sensed bya point sink positioned in a being at rest point y at the given time tby applying Equation (36) given below:

$\begin{matrix}{{{J_{y^{\prime}\rightarrow y}^{\Psi}\left( {t,y} \right)} = {{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{2{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}}{\Psi_{i\; n}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{dy}^{\prime}}};} & (36)\end{matrix}$

and defining the rate of the property vector flux, J_(FS→y) ^(Ψ), inpoint y at the given time t, which is originated from precedingcollisions within entire space occupied by the model gas, includesapplying Equation (37) given below, which is obtained by integratingEquation (36) over the volume of the model gas system:

$\begin{matrix}{{J_{{FS}\rightarrow y}^{\Psi} = {\frac{1}{2}{\int_{y_{b\; 2}}^{y_{b\; 1}}{{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}{\Psi_{i\; n}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{dy}^{\prime}}}}},} & (37)\end{matrix}$

where Ψ_(in)(t_(i)′,y′,t,y) is the property carried by each of theconverging ballistic particles at the moment of the entering the volumesurrounding point y at the given time t.

FIG. 6 shows the model gas system schematically with a uniform model gasflow along the x-direction, in which a main control volume 614 of volumedV=Δxdy is at y within the model gas system confined by two parallelplates at y_(b1) and y_(b2). In the example of FIG. 6, the acceleration,600, is directed in the negative to y-axis direction. The arrow 601 or603 indicates schematically flux expressed by Equation (36) andassociated with divergence of particles from an auxiliary control volume602 or from an auxiliary control volume 604, respectively.

The step of defining the net rate of property influx per unit volume,B_(in) ^(Ψ0_F), formed by the flow of ballistic particles carryingproperty and converging from the space occupied by the model gas in thegeneral non-moving point of the space occupied by the model gas at thegiven time further includes a step of representing B_(in) ^(Ψ0_F) byapplying Equation (38) given below, which is obtained by executing

$\frac{\partial}{\partial y}$

operator to the one-dimensional vector field of Equation (37):

$\begin{matrix}{{B_{i\; n}^{{\Psi\_}F}\left( {y,t} \right)} = {{- \frac{1}{2}}\frac{\partial}{\partial y}{\int_{y_{b\; 2}}^{y_{b\; 1}}{{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}{\Psi_{i\; n}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{{dy}^{\prime}.}}}}} & (38)\end{matrix}$

1.2.2 A Combined Analytical Tool for Computing the Net Rate of PropertyInflux in the General Non-Moving Point at the Given Time, which isFormed by Converging Ballistic Particles Originated from Gas-SolidInterfaces

FIG. 7 is a schematic view of the model gas system confined betweenparallel plates, shown to illustrate a variety of ballistic trajectoriesof particles in the model gas system having gas-solid interfaces withmixed diffuse and specular components of particle scattering. In theexample of FIG. 7, the one-dimensional model gas system is confinedbetween spaced a top plate 708 and a bottom plate 709. A channel formedby spaced plates has an inlet opening on the left side, shown in theschematic by line AB, and the outlet opening on the right side, shown byline CD. The schematic shows a variety of types of ballistictrajectories in a system with gas-solid interfaces revealing bothdiffuse and specular scattering components, which are arranged forillustration and do not intend to limit the invention. The first type,ballistic trajectory 701 shows a way for transporting property carriedby the particle leaving an auxiliary control volume 706 in y′ and comingdirectly in point y of a main control volume 707. The second type,ballistic trajectories 702 or 703 show ways for transporting propertycarried by particles leaving after diffuse scattering from plate 708 andplate 709, respectively, and coming directly in point y of the controlvolume. The third type, ballistic trajectory 704 with the final segment704A shows a possible way for transporting a property carried by theparticle leaving auxiliary control volume 706 toward plate 708, then isspecularly scattered back into model gas toward plate 709 and may followone or more of specular scatterings from gas-solid interfaces until itfinally makes final ending collision in point y of main control volume707. The fourth type, ballistic trajectory 705 and the final trajectorysegment 705A shows a way for transporting a property carried byparticles leaving by diffusive scattering plate 709 toward plate 708,then is specularly scattered back into model gas toward plate 709 andmay follow one or more of specular scatterings from gas-solid interfacesuntil it finally makes following final collision in point y being withinthe main control volume. Specifically, according to the example of FIG.7, determining total influx of property delivered into main controlvolume 707 by traveling particles includes the integration of propertyintake according to 701, 702, 703, 704A and 705A ending segments ofballistic trajectories.

It is well known that since the solid surface is not ideal, it containsboth diffuse and specular components of particle scattering. The diffusescattering implies that each particle that hits a plate transfers itsmomentum/energy or other property completely to the plate and that allparticles diffusively scattered from the plate leave it acceptingstreamwise velocity and temperature or other property of the plate. Ifparticle collisions are specular the incident on a plate property, whichis carried by the particle such as energy, will be returned entirelyinto the model gas system.

By the method, the diffuse reflection is characterized by the propertyaccommodation coefficient 0<σ≤1, and the specular reflection ischaracterized by 1−σ. Each particle diffusively scattered from the plate1 leaves it accepting one or more properties intrinsic to the gas-solidinterface of the plate such as the thermal velocity.

For clarity, all gas-solid interfaces bounding the model gas flow areconsidered to have identical gas-solid interfaces, which arecharacterized by unchangeable property accommodation coefficient, σ. Itis to be understood this limitation is made without departing from theclaiming subject.

For diffuse scattering with no adsorption effects after colliding withthe surface, the particle flux per unit time diffusively scattered backfrom the gas-solid interface into the model gas system equals to theparticle flux per unit time incident on the surface of the boundary witha negative sign. Here we have also accepted that particles diffusivelyscattered back from the gas-solid interface will accommodate a thermalvelocity and other properties at the moment of the contact with thesurface of the boundary at a particular place of the boundary. Also, wehave recognized that only a half of the particles near scatteringsurface had a recent departure from the surface with the magnitude ofvelocity equal to the magnitude of thermal velocity corresponding to thetemperature of the reflecting gas-solid interface.

The step S103 of defining the net rate of property influx per unitvolume from the model gas system further includes a step of defining anet rate of total property influx per unit volume from a gas-solidinterface exhibiting mixed diffuse and specular particle scattering,

where the net rate of total property influx per unit volume from thegas-solid interface is formed by way of a second plurality of convergingballistic particles, each of the second plurality of convergingballistic particles is selected from the plurality of convergingballistic particles by the ballistic trajectory having the startingpoint in one of a plurality of points of original collisions on thegas-solid interface, and

where each of the plurality of points of original collisions resulted indiffuse particle scattering from the gas-solid interface is treated as aheterogeneous point source.

FIG. 8 is a block diagram of an alternative approach employed to definemajor steps needed to calculate the net rate of property influx per unitvolume from the diffuse gas-solid interfaces of the boundary in thegeneral non-moving point y at the given time t.

In some embodiments, the step S103B in the example of FIG. 1A ofdefining the net rate of property influx per unit volume in the generalnon-moving point at the given time from the gas-solid interface of theboundary includes steps of:

identifying the ballistic trajectory, for each of the second pluralityof converging ballistic particles (S103B1),

where the velocity of the heterogeneous point source for each of thesecond plurality of converging ballistic particles is assigned to beequal to the velocity of the gas-solid interface in a correspondingpoint of original collisions on the gas-solid interface at the time ofdiffuse particle scattering;

defining the probability of traveling along the ballistic trajectory(S103B2);

defining a rate of collisions per unit area in one of a plurality ofheterogeneous points on the gas-solid interface at the time of diffuseparticle scattering (S103B3),

where the ballistic trajectory identifies the time of diffuse particlescattering for each of the second plurality of converging ballisticparticles;

defining a net rate of property flux per unit area in the general pointat the given time from one of a plurality of heterogeneous point sources(S103B4);

defining a net rate of total property flux in the general point at thegiven time from the plurality of heterogeneous point sources on one ormore of gas-solid interfaces (S103B5); and

defining the net rate of total property influx per unit volume in thegeneral point at the given time by applying the divergence operator tothe net rate of total property flux (S103B1).

In the example of FIG. 6, plate 615 and plate 616 confine across y-axisa model gas flow along the x-axis and surfaces 611 and 612 positioned ata distance Δx bound a portion of the model gas system along the x-axis.In the model gas system having gas-solid interfaces with purely diffusescattering of particles, each particle originated from preceding diffusescatterings from the interfaces delivers in the main control volume 614some property obtained from the location and time of the originaldiffuse scattering (ballistic trajectory 606 from plate 616 andballistic trajectory 605 from plate 615 of FIG. 6).

To calculate property fluxes in a originated from diffuse particlescattering from the gas-solid interfaces of two stationary parallelplate s it is needed to determine the departing time t_(ib1)′ andt_(ib2)′ from plate 615 and plate 616, respectively, which provide anopportunity for departing the plate particles to arrive in point y attime t. We have also recognized that departing time from plate 615 orplate 616 is needed to indicate the time at which the departing particlehas obtained a particular transported property. Referring again to theexample of FIG. 6, determining the departing times t_(ib1)′ and t_(ib2)′from plate 615 and plate 616, respectively, includes applying formulas:

t _(ib1) ′=t−φ _(ib1)′ and t _(ib2) =t−φ _(ib2),  (39)

where φ_(ib1)′ and Φ_(ib2)′ are ballistic traveling times from plate 615and plate 616, respectively, to point y, which can be determined bysolving corresponding equations of projectile motion with respect tocorresponding traveling times φ_(ib1)′ and φ_(ib2)′.

Referring again to FIG. 6 showing schematically model gas flow confinedby two being at rest parallel plates plate 615 and plate 616 at y_(ib1)′and φ_(ib2)′ and having normal vectors n_(b1) and n_(b2), respectively,upon recognizing that the property flux in a general non-moving point yat time t, is originated from the original collisions with gas-solidinterfaces of plate 615 and plate 616, the property vector fluxes perunit time from plate 615 along ballistic trajectory 605 and plate 616along ballistic trajectory 606 is calculated from Equations (40) and(41), respectively, given below:

$\begin{matrix}{{J_{y_{b\; 1}\rightarrow y}^{\Psi}\left( {t,y} \right)} = {\sigma \; {Q_{{ib}\; 1}\left( {t,t_{{ib}\; 1}^{\prime}} \right)}{Z_{b\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}\frac{v_{b\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1},t,y} \right)}{v_{{Tb}\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}{\Psi_{\lambda \; b\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1},\lambda,\phi_{{ib}\; 1}} \right)}}} & (40) \\{\mspace{79mu} {and}} & \; \\{{{J_{y_{vb2}\rightarrow y}^{\Psi}\left( {t,y} \right)} = {\sigma \; {Q_{{ib}\; 2}\left( {t,t_{{ib}\; 2}^{\prime}} \right)}{Z_{b\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2}} \right)}\frac{v_{b\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2},t,y} \right)}{v_{{Tb}\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2}} \right)}{\Psi_{\lambda \; b\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2},\lambda,\phi_{{ib}\; 2}} \right)}}},} & (41)\end{matrix}$

where Z_(b1)(t_(ib1)′,y_(b1)) or Z_(b2)(t_(ib2)′,y_(b2)) is the totalnumber of collisions per unit area per unit time on surfaces positionedin y_(b1) or y_(b2) at time t_(ib)′ or t_(ib2)′, respectively;Q_(ib1)(t,t_(ib1)′) or Q_(ib2)(t,t_(ib2)′) is the probability ofballistic traveling along the trajectory beginning from y_(b1) or y_(b2)at time t_(ib1)′ or t_(ib2)′, respectively, allowing each of thetraveling particles to arrive in point y at time t,Ψ_(λb1)(t_(ib1)′,y_(b1),λ,φ_(ib1)) or Ψ_(λb2)(t_(ib2)′,y_(b2),λ,φ_(ib2))is property obtained by the traveling particle from plate 615 or plate616, respectively; v_(b1)(t_(ib1)′,y_(b1),t,y) orv_(b2)(t_(ib2)′,y_(b2),t,y) is the velocity of the traveling particletargeting point y at the time t which is originated from the particle'soriginal diffuse scattering from point y_(b1) or y_(b2) at time t_(ib1)′or t_(ib2)′, respectively; v_(Tb1)(t_(ib1)′,y_(b1)) orv_(Tb2)(t_(ib2)′,y_(ib2)) is the magnitude of the thermal velocity ofthe particle leaving the gas-solid interface of plate 615 or plate 616at time t_(ib1)′ or t_(ib2)′, respectively.

In the model gas confined between two infinite parallel plates, theparticles diffusively scattered from one plate have a chance to bedelivered to the other plate without intermediate collisions in themodel gas space.

In some embodiments, referring again to FIG. 6, defining the totalnumber of collisions per unit area per unit time on a gas-solidinterface of one of two stationary parallel plates confining the modelgas flow and revealing the purely diffuse mechanism of the particlescatterings, i.e., σ=1, includes steps of:

(1) determining particle vector flux, N₂ ¹, 607 incident on plate 615 attime t_(ib1)′, which is originated from diffuse scattering from thegas-surface interface of plate 616, by calculating Equation (42) givenbelow, which is obtained by assigning Ψ_(b1)=1, y=y_(b1), u=0, andt=t_(ib1)′ in Equation (41):

$\begin{matrix}{{N_{2}^{1} = {{Q_{ib}\left( {t_{{ib}\; 1}^{\prime},{t_{{ib}\; 1}^{\prime} - \phi_{{b\; 2} - {b\; 1}}}} \right)}{Z_{b}\left( {{t_{{ib}\; 1}^{\prime} - \phi_{{b\; 2} - {b\; 1}}},y_{b\; 2}} \right)}\frac{v_{b\;}\left( {{t_{{ib}\; 1}^{\prime} - \phi_{{b\; 2} - {b\; 1}}},y_{b\; 2},t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}{v_{{Tb}\; 1}\left( {{t_{{ib}\; 1}^{\prime} - \phi_{{b\; 2} - {b\; 1}}},y_{b\; 2}} \right)}}},} & (42)\end{matrix}$

where φ_(b2−b1) is traveling time from plate 616 to plate 615;

(2) determining the particle vector flux per unit time, N₂₁ ¹, shownschematically as 608 incident on plate 615 at time t_(ib1)′, which isassociated with the incident on the gas-solid interface of plate 615 ofballistic particles from the gas in the range from y_(b2) to y_(b1), iscalculated from Equation (43) given below, which is obtained byassigning t=t_(ib1)′, t_(i)′=t_(ib1)′, −φ_(i1), y=y_(b1), Ψ_(in)=1, andu=0 in Equation (37):

$\begin{matrix}{{N_{21}^{1} = {\frac{1}{2}{\int_{y_{b\; 2}}^{y_{b\; 1}}{{Q_{i}\left( {t_{{ib}\; 1}^{\prime},{t_{{ib}\; 1}^{\prime} - \phi_{i\; 1}^{\prime}}} \right)}{Z_{V}\left( {{t_{{ib}\; 1}^{\prime} - \phi_{i\; 1}^{\prime}},y^{\prime}} \right)}\frac{v\left( {{t_{{ib}\; 1}^{\prime} - \phi_{i\; 1}^{\prime}},y^{\prime},t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}{v_{T}\left( {{t_{{ib}\; 1}^{\prime} - \phi_{i\; 1}^{\prime}},y^{\prime}} \right)}{dy}^{\prime}}}}},} & (43)\end{matrix}$

where φ_(i1)′ is traveling time between y′ and plate 615,Q_(i)(t_(ib1)′,t_(ib1)′−φ_(i1)′) is the survival probability of theparticle ballistic traveling along the ballistic trajectory from thepoint of the preceding collision on the boundary at y_(b2) to theboundary at y_(b1) at time t_(ib1), which, is expressed as follows

Q _(i)(t _(ib1) ′,t _(ib1)′−φ_(i1)′)exp(−∫_(t) _(ib1) _(′) _(−φ) _(i1)_(′) ^(t) ^(ib1) ^(′) P _(c)(y({tilde over (t)},t _(ib1)′−φ_(i1)))v_(rel) y({tilde over (t)},t _(ib1)′−φ_(i1)))d{tilde over (t)});  (44)

(3) determining the impingement rate, J_(FS) ^(ib1), or the totalparticle vector flux per unit time incident at time t_(ib1)′ on thesurface of the plate 615 at y_(b1), by calculating Equation (45) givenbelow, which is obtained by summing Equations (42) and (43):

J _(FS) ^(yb1) =N ₂₁ ¹ +N ₂ ¹;  (45)

(4) determining the rate of collisions Z_(b1) at time t_(ib1)′ on asurface of plate 615 having normal n_(b1) or the impingement rate iscalculated from Equation (46) given below:

Z _(b1) =J _(FS) ^(yb1) n _(b1)  (46)

Analogously, the impingement rate, J_(FS) ^(yb2), or the total vectorflux of particles per unit area incident at time t_(ib2)′ on the surfaceof the plate at y_(b2) is calculated as follows

J _(FS) ^(yb2) =N ₁₂ ² +N ₁ ²,  (47)

where N₁ ² is particle flux 609 on plate 616 at time t_(ib2)′, which isformed by ballistic traveling particles diffusively scattered from theplate 615, N₁₂ ² is the particle flux 610 incident at time t_(ib2)′ onthe surface of plate 616;

Finally, the rate of collisions Z_(b2) at time t_(ib2)′ on a surface ofplate 616 having normal n_(b2) or the impingement rate is calculatedfrom Equation (48) given below:

Z _(b2) =−J _(FS) ^(yb2) n _(b2).  (48)

The step S103B of defining the net rate of property influx per unitvolume from the diffuse gas-solid interface of the boundary in the samegeneral non-moving point at the same time includes, for example, a stepof quantifying property supplied by converging ballistic particlesoriginated from diffuse scatterings from gas-solid interfaces of plate615 having normal n_(b1)=−1 and positioned at y_(b1) and plate 616having normal n_(b2)=1 and positioned at y_(b2) by applying Equations(49) and (50), respectively, given below, which is obtained by executing

$\frac{\partial}{\partial y}$

operator to the one-dimensional vector field of Equations (40) and (41)and multiplying by a negative sign, respectively:

$\begin{matrix}{B_{i\; n}^{\Psi \; {{\lambda\_}{bd}}\; 1} = {- {\frac{\partial}{\partial y}\left\lbrack J_{y_{b\; 1}\rightarrow y}^{\Psi} \right\rbrack}}} & (49) \\{and} & \; \\{{B_{i\; n}^{\Psi \; {{\lambda\_}{bd}}\; 2} = {- {\frac{\partial}{\partial y}\left\lbrack J_{y_{b\; 2}\rightarrow y}^{\Psi} \right\rbrack}}},} & (50)\end{matrix}$

where B_(in) ^(Ψλ_bd1) and B_(in) ^(Ψλ_bd2) are the net rate of propertyinflux per unit volume from the diffuse gas-solid interface of theboundary formed by converging ballistic particles originated fromdiffuse scatterings from gas-solid interfaces of plate 615 and plate616, respectively.

For the model gas system established to be a bounded model gas systemproviding a model gas flow passage through a channel formed by aboundary confining the channel and being in contact with the model gasthrough a model gas-solid interface exhibiting mixed diffuse andspecular particle scattering, the step S103 in the example of FIG. 1A ofdefining the net rate of property influx per unit volume in the generalnon-moving point of the space occupied by the model gas at the giventime from the surrounding model gas system in some embodiments includessteps of:

defining the net rate of property influx per unit volume in the generalnon-moving point at the given time from the surrounding model gas, whichis formed by the converging ballistic traveling particles originatedfrom the original collisions within the space occupied by the model gas(S103A in the example of FIG. 3);

defining the net rate of property influx per unit volume in the generalnon-moving point at the given time from the diffuse gas-solid interfaceof the boundary, which is formed by the converging ballistic travelingparticles originated from the original scatterings from the gas-solidinterface exhibiting diffuse particle scattering (S103B);

defining the net rate of property influx per unit volume in the generalnon-moving point at the given time via the mixed diffuse-speculargas-solid interfaces of the boundary (S103C),

where the value of property delivered in the general non-moving point atthe time by each of the converging ballistic traveling particles in themodel gas system is evaluated regarding whether the value of theproperty is conserved or changed because of aging such as change of masswith internal radioactive decay during the ballistic traveling time,

where the value of property delivered in the general non-moving point atthe given time by each of the converging ballistic particles in themodel gas system is evaluated regarding whether the value of theproperty is modified because of interaction with an external field suchas change of the velocity with a gravitational field during theballistic traveling time;

defining the net rate of property generation per unit volume in thegeneral non-moving point originated from the effect of surrounding suchas the model gas pressure gradient or the temperature gradient (S103D);and

summing influxes of steps S103A, S103B, 103C, and S103D according tostep S103E.

The step S103D of defining the net rate of property generation per unitvolume in the general non-moving point originated from the effect ofsurrounding such as the model gas pressure gradient or the temperaturegradient includes, for example, a step of adding contribution oftransported quantity within the model gas system, which is resulted fromtransmitting property across interfaces separating the model gas of thesystem and the model gas of the surroundings such as the surfaces ofinlet and outlet openings.

In some embodiments of a bounded model gas system having an open-channelgeometry providing fluid flow passage through a channel having an inletopening and an outlet opening, which is defined by approximating to anidentical bounded model gas system providing model gas flow passagethrough a channel formed by a boundary confining the channel and beingin contact with the model gas through a gas-solid interface, additionalgeneration of transported quantity within the model gas system isresulted from transmitting property across the surfaces of inlet andoutlet openings. Symbolically, the net rates of property generationoriginated from the effect of the surrounding (model gas pressuregradient, temperature gradient) is expressed as Q_(T).

If the property transport is the momentum transport, then the generationterm in a confined model gas system results from forces acting on themodel gas in the form of the external pressure force, P, applied acrossthe surface bounding the system. We have also recognized that normalgradient pressure force acting on the model gas and prompting model gasflow represents the effect of surrounding on the momentum change.Particularly, referring to FIG. 6, if external pressure P_(x) is appliedto surface 611 in the positive x-direction and external pressureP_(x+Δx) is applied to surface 612 in the negative x-direction and thelength of the model gas system is Δx, then the negative pressuregradient along the x-axis,

${- \frac{P_{x + {\Delta \; x}} - P_{x}}{\Delta \; x}},$

is treated as the force per unit volume motivating the model gasmovement toward the resulting force, i.e.:

$\begin{matrix}{{\overset{.}{Q}}_{\Psi} = {- {{\lim \left( \frac{P_{x + {\Delta \; x}} - P_{x}}{\Delta \; x} \right)}_{{\Delta \; x}\rightarrow 0}.}}} & (51)\end{matrix}$

If property transport is energy transport, the generation terms resultfrom energy flow through the main control volume because of inducedtemperature gradients.

The embodiments of the following part of the specification are explainedwhere property acquired by the particle in the point of the precedingcollision or preceding diffuse scattering is not aged during aparticle's ballistic traveling time, which implies that transportedproperty is conserved during ballistic traveling and aging coefficientλ=0

FIG. 9 is a block diagram of an alternative approach employed to definemajor steps of step S103C needed to calculate the net rate of propertyinflux per unit volume via one or more of mixed diffuse-speculargas-solid interfaces in the general non-moving point y at the given timet. The step S103C includes steps of:

defining the net rate of property influx per unit volume formed by theconverging ballistic particles providing delivery of property obtainedby each from a corresponding to it point source in the space occupied bythe model gas through at least one preceding sequential specularscattering from the gas-solid interface (S103C1);

defining the net rate of property influx per unit volume formed byconverging ballistic particles delivering property obtained by each froma corresponding to it point source on the gas-solid interface through atleast one preceding sequential specular scattering from said gas-solidinterface (S103C2),

where point source strength on the gas-solid interface having bothdiffuse and specular particle scattering components must be proportionalto the property accommodation coefficient σ;

summing influxes of steps S103C1 and S103C2 according to step S103C3.

Further, the embodiments of this and the following sub-sections dealingwith mixed diffuse and specular scattering are explained in a case ofone-dimensional incompressible steady-state model gas flow between beingat rest infinite parallel plates having gas-solid interfaces with mixedboth diffuse and specular particle scattering and at the uniformtemperature implying that the thermal velocity of the particles is aconstant. The present invention, however, is not limited to the case butcan be preferably applied to numerous variations and modifications.

FIG. 10 is a schematic view of the one-dimensional model gas systembetween parallel infinite plates, shown to illustrate a schematic ofdefining diffuse and specular components of the impingement rate on thegas-solid interfaces the plates. Specifically, the impingement rate onthe gas-solid interface of plate 1009 is associated with the particlesincident on plate 1009, which are originated from preceding collisionswithin of both the model gas confined between the interfaces of twoparallel plates, plate 1009 at y_(b1) and plate 1010 at y_(b2) with theproperty accommodation coefficient on both gas-solid interfaces, σ. Theschematic shows also various particle paths from the model gas system onthe plates and paths between the confining plates. In the example ofFIG. 10, which shows the model gas system having gas-solid interfaceswith mixed diffuse and specular components of particle scattering, avariety of particle paths are involved in the transport of the particleson a gas-solid interface, with each path being characterized by theprobability of the ballistic traveling. Defining the impingement rate ona specific location of the gas-solid interface at the given timeincludes, for example, steps of:

defining a bulk component of an impingement rate formed by particleschosen from the plurality of converging ballistic particles originatedfrom the original collisions in space occupied by the model gas(trajectory 1001);

defining a bulk-specular component of the impingement rate formed byparticles chosen from the plurality of converging ballistic particlesoriginated from the original collisions in space occupied by the modelgas (trajectory 1003) and having at least the last preceding specularscattering (trajectory 1005);

defining a diffuse component of the impingement rate formed by particleschosen from the plurality of converging ballistic particles originatedfrom one or more of heterogeneous point sources (trajectory 1002); and

defining a diffuse-specular component of the impingement rate formed byparticles chosen from the plurality of converging ballistic particlesoriginated from original diffuse scatterings from one or more ofheterogeneous point sources (trajectory 1004) and having at least thelast preceding specular scattering (trajectory 1006); and

summing the bulk component, the bulk-specular component, the diffusecomponent, and the diffuse-specular component of the impingement rate,

where each of a plurality of impinging particles is selected by animpinging trajectory allowing to target the gas-solid interface in aspecific location of the gas-solid interface at a specific time.

In the model gas confined between two infinite parallel plates, some ofthe particles having original collisions within the model gas volume andsome particles diffusively scattered from plate 1010 have a chance totravel to the other plate 1009 without intermediate collisions in themodel gas space (according to ballistic trajectories 1001 and 1002,respectively). Some of particles having original collisions within themodel gas volume and some particles diffusively scattered from plate1009 have a chance to be delivered to plate 1010 without intermediatecollisions in the model gas space (according to ballistic trajectories1003 and 1004, respectively) and then specularly scattered again back toplate 1009 (according to ballistic trajectories 1005 and 1006,respectively. Also, since both gas-solid interfaces confining the modelgas have specular components, some particles scattered back have achance to bounce between the plates (as shown in FIG. 10 schematicallyby bouncing cycles 1007 and 1008.

Since the model gas system is kept at the uniform constant temperature,all particle fluxes per unit time within the model gas system arestable. Under the method, property accommodation coefficient σ on thegas-solid interface having both diffuse and specular particle scatteringis treated as the probability of incident particles to be “absorbed”into the surface with further its “evaporation” back with propertiesbeing acquired by the particle at the gas-solid interface at the momentof “absorption”.

Referring again to FIG. 10, upon recognizing that the total number ofparticles colliding with a surface of a plate per unit area per unittime or the impingement rate is formed by the impinging particlesaccording to ballistic trajectories described above and bouncingballistic trajectories (shown as the first and the second bouncingcycles 1007 and 1008, respectively), the impingement rates or theparticle vector fluxes on plate 1009, J_(FS) ^(yb1), and plate 1010,J_(FS) ^(yb2) are obtained from solving a system of Equations (52) and(53), given below

J _(FS) ^(yb1)=[N ₂₁ ¹−(1−σ)N ₁₂ ² exp(−P _(c) H)−σJ _(FS) ^(yb2) exp(−P_(c) H)+σ(1−σ)J _(FS) ^(yb1) exp(−2P _(c) H)]

  (52)

and

J _(FS) ^(yb2)=[N ₁₂ ²−(1−σ)N ₂₁ ¹ exp(−P _(c) H)−σJ _(FS) ^(yb1) exp(−P_(c) H)+σ(1−σ)exp(−2P _(c) H)J _(FS) ^(yb2)]

  (53)

where, in the equations above, the first terms are N₂₁ ¹ and N₁₂ ²,which are particle fluxes per unit time incident on plate 1009 and plate1010 formed by the particles following ballistic trajectories 1001 and1003, respectively, the second terms are −(1−σ)N₁₂ ² exp(−P_(c)H) and−(1−σ)N₂₁ ¹ exp(−P_(c)H), which are specular scattered back fluxes perunit time from plate 1010 and incident on plate 1009 (trajectory 1005)and from plate 1009 (top) and incident on plate 1010, respectively, thethird terms are −σJ_(FS) ^(yb2) exp(−P_(c)H) and −σJ_(FS) ^(yb1)exp(−P_(c)H), which are the particle fluxes per unit time emitteddiffusively from point sources, point source strength of each isdirectly proportional to σ, of plate 1010 and incident on plate 1009)(trajectory 1002) from plate 1009 and incident on plate 1010,respectively, and the fourth terms are σ(1−σ)J_(FS) ^(yb1) exp(−2P_(c)H)and +σ(1−σ)exp(−2P_(c)H)J_(FS) ^(yb2), which are the particle fluxes perunit time scattered diffusively from plate 1009 to plate 1010 and thenspecularly scattered back to plate 1009 (trajectory 1006) and from plate1010 to plate 1009 and then specularly scattered back to plate 1010,respectively. In the equations above,

is the bouncing coefficient or a coefficient characterizing the totalimpact of the sequential specular bouncing on the impingement rates onboth plate 1009 and plate 1010. Under the method, the bouncingcoefficient is calculated from Equations (54) given below:

$\begin{matrix}{{{\mathbb{N}} = {{\sum\limits_{n = 0}^{\infty}{\left( {1 - \sigma} \right)^{2n}{\exp \left( {{- 2}{nP}_{c}H} \right)}}} = \frac{1}{1 - {\left( {1 - \sigma} \right)^{2}{\exp \left( {{- 2}P_{c}H} \right)}}}}},} & (54)\end{matrix}$

where (1−σ)² exp(−2P_(c)H) is the impact of one cycle of bouncingstarting from both plate 1009 or plate 1010.

Referring again to FIG. 10, upon recognizing that, the model gas isconfined between two infinite parallel plates having identicalproperties on the gas-solid interfaces, so that

N ₁₂ ² =−N ₂₁ ¹ and J _(FS) ^(yb2) =−J _(FS) ^(yb1),  (55)

and N₂₁ ¹ is calculated according to Equation (56) given below, which isobtained by integrating Equation (43) with P_(c)=constant andZ_(V)=constant at the uniform temperature:

$\begin{matrix}{{N_{21}^{1} = {\frac{Z_{V}}{2P_{c}}\left\lbrack {1 - {\exp \left( {{- P_{c}}H} \right)}} \right\rbrack}},} & (56)\end{matrix}$

then the impingement rate on plate 1009 or the rate of collisions,J_(FS) ^(yb1) is calculated from Equations (57), given below, which isobtained by substitution of Equations (56) and (55) in Equation (52) or(53) and executing algebra operations:

$\begin{matrix}{Z_{b} = {J_{FS}^{{yb}\; 1} = {{- J_{FS}^{{yb}\; 2}} = {+ {\frac{Z_{V}}{2P_{c}}.}}}}} & (57)\end{matrix}$

FIG. 11 is a schematic view of the one-dimensional model gas system,shown to illustrate a method for analytical representation of the netrate of property influx in a general non-moving point y, which is formedby converging ballistic particles, each having at least one specularscattering from a gas-solid interface of plate 1105 or plate 1106.

In a stable or steady-state model gas system, any propertycharacterizing model gas flow depends on only the position in space, y′,but not time, i.e., Ψ(y′) and properties characterizing the gas-solidinterfaces of plate 1105 positioned at y_(b1)=H and plate 1106positioned at y_(b2)=0, Ψ_(b1) and Ψ_(b2), respectively, are constants.The embodiments of this part of the specification are explained wherethe values of the property on the plate is assigned to be zero, namely,Ψ_(b1)=Ψ_(b2)=0. This will cause eliminating fluxes associated withconverging ballistic particles chosen from the group of particlesoriginated from original diffuse scatterings from plate 1105 and plate1106. This limitation is made for clarity without limiting the presentinvention.

In the stable system with the incompressible model gas flow, theprobability of ballistic traveling, Q_(i), between points y′ and y iscalculated from Equation (58) given below, which is obtained bysubstitution of v_(rel)=v_(T), {tilde over (y)} of Equation (10) withg=0 and u=0, and P_(c)=constant in Equation (21):

Q _(i)(y,y′)=exp(−P _(c) |y−y′|).  (58)

Concerning the example of FIG. 11, the net rate of property influx,B_(in) ^(Ψ0_F_bs), formed by converging ballistic particles chosen fromthe group of particles originated from initial collisions within themodel gas volume and having at least the last preceding specularscattering from gas-solid interfaces of both plate 1105 and plate 1106is calculated from Equation (59) given below

B _(in) ^(Ψ0_F_bs) =B _(in) ^(Ψ0_F_bs1) +B _(in) ^(Ψ0_F_bs2),  (59)

where B_(in) ^(Ψ_F_bs1) and B_(in) ^(Ψ_F_bs2) are the net rate ofproperty influxes associated with the last specular scatterings fromplate 1105 and plate 1106, respectively. Specifically, B_(in) ^(Ψ_F_bs1)is formed by particles having ballistic trajectories 1101 and 1102,which are continued by segmented trajectory 1103 with the lasttrajectory segment 1104 released from plate 1105, as shown in FIG. 11and is calculated from Equation (60) given below:

$\begin{matrix}{{B_{i\; n}^{\Psi \; 0{\_ F}{\_ bs}\; 1} = {\frac{\partial}{\partial y}\left\{ {\left\lbrack {{\left\lbrack {1 - \sigma} \right\rbrack J_{21}^{1}} - {\left\lbrack {1 - \sigma} \right\rbrack^{2}J_{12}^{2}{\exp \left( {{- P_{c}}H} \right)}}} \right\rbrack {\mathbb{N}}\; {\exp \left( {- {P_{c}\left( {H - y} \right)}} \right)}} \right\}}},} & (60)\end{matrix}$

where J₂₁ ¹ and J₁₂ ² are the property fluxes per unit time, which areformed by ballistic particles incident on plate 1105 and plate 1106,respectively, and

is calculated from Equation (54).

The property vector fluxes J₂₁ ¹ and J₁₂ ² are calculated from Equations(61) and (62) given below:

J ₂₁ ¹=½z _(V)∫₀ ^(H)Ψ(y′)exp(−P _(c) y′)dy′  (61)

and

J ₁₂ ²=−½ exp(−P _(c) H)z _(V)∫₀ ^(H)Ψ(y′)exp(P _(c) y′)dy′,  (62)

respectively. Analogously, using Equation (60) given above as atemplate, we obtain

$\begin{matrix}{B_{i\; n}^{\Psi \; 0{\_ F}{\_ bs}\; 2} = {\frac{\partial}{\partial y}{\left\{ {\left\lbrack {{{- \left\lbrack {1 - \sigma} \right\rbrack}J_{12}^{2}} + {\left\lbrack {1 - \sigma} \right\rbrack^{2}J_{21}^{1}{\exp \left( {{- P_{c}}H} \right)}}} \right\rbrack {\mathbb{N}}\; {\exp \left( {{- P_{c}}y} \right)}} \right\}.}}} & (63)\end{matrix}$

Next, to define and formulate the net rate of property efflux per unitvolume into surroundings from the general non-moving point y at thegiven time, t, the linear dimensions of the main control volumesurrounding point y are selected to be sufficiently small for avoidingtwo and more consecutive collisions of the same particle within the maincontrol volume.

1.2.3 A Combined Analytical Tool for Computing the Net Rate of PropertyEfflux from the General Non-Moving Point at the Given Time, which isFormed by Diverging Ballistic Particles

FIG. 12 is a schematic shown to illustrate a method for analyticalrepresentation of the net rate of property efflux from a stationarypoint in one-dimensional space occupied by the model gas. According tothe schematic, a particle 1200 in point y at the given time t divergesfrom the point source in this point and achieves velocity v₊ in point y′at time t_(a)′ while crossing a control surface 1201 or 1202 of acontrol volume 1203 surrounding the point source. We assume thatproperty, Ψ_(λ+) carried by the diverging ballistic particle at timet_(a)′ of crossing in point y′ enclosing control surface 1201 or 1202was initially obtained by the particle at the moment of the originalcollision in point y at time t. Specifically, in a point source in y attime t, (1) each particle diverging from a general non-moving point y atthe time, t, has acquired the thermal velocity v_(T)(t,y), mass flowvelocity u(t,y), and one or more specific properties, Ψ_(λ+)(t,y), ofparticular values including mass, momentum, and energy; (2) eachtraveling particle is being simultaneously in the states of bothdivergence from the point of preceding collision and convergence in thepoint of sequential collision; (3) equal number of converging anddiverging ballistic particles is maintained at any point of collisionsat any time, which provides maintaining continuity in the model gas. Inthe example of FIG. 12, the observer's coordinate system is orientedwith y-axis directed along the negative direction of the appliedacceleration field 1204 with acceleration g for each particle. Theparticle having a collision in point y at time t acquires a magnitude ofthermal velocity, v_(T)(t,y) and mass flow velocity u(t,y), which areintrinsic to the nearby volume of the model gas surrounding point y attime t. After a collision in point y at time t, the particle may move inany direction because of the arbitrary nature of the unit vector n_(+c),thus forming an expansion zone 1203.

The step S104 of defining the net rate of property efflux per unitvolume from the general point of the plurality of non-moving points atthe given time in the model gas system in one-dimensional configurationincludes steps of:

for each of the plurality of diverging particles, identifying theballistic trajectory starting from the general point of the plurality ofnon-moving points at the given time and ending in one of the pluralityof points in space surrounding the general point of the plurality ofnon-moving points;

defining the probability of traveling along the ballistic trajectoryfrom the point of the divergence y to point y′ on the control surface1201 or 1202;

defining a vector field of a property flux per unit area in the endingpoint of one of the plurality of points in space surrounding the generalpoint;

defining the net rate of property efflux per unit volume by applying adivergence operator to the vector field of the property flux per unitarea in the ending point; and the shrinking volume of space surroundingthe general point to the general point.

The step of identifying a trajectory and trajectory characteristics, foreach particle diverging from the general non-moving point at the time,includes, for example, steps of:

defining a diverging ballistic trajectory of each particle divergingfrom the general non-moving point y at time, t;

determining traveling time needed, for each particle diverging fromnon-moving point y at the given time, t, to cross a control surfaceenclosing point y in point y′ of control surface;

defining an appropriate instant unit vector directing thermal velocitycomponent along the diverging ballistic trajectory of each particlediverging from the general non-moving point y at time t;

defining the velocity vector of each particle at the moment of crossingthe control surface in point y′.

1.2.3.1 Analytical Tool for Computing the Trajectory and Velocity of aBallistic Movement of the Particle Diverging from the General Non-Movingat the Given Time

The step of defining the diverging ballistic trajectory of each particlediverging from the general non-moving point y at time t includes a stepof formulating position vector {tilde over (y)}′ 1205 or 1206 of theparticle on trajectory 1207 or 1208, respectively, and velocity vector{tilde over (v)}′ (not shown) at time {tilde over (t)} by applyingEquations (64) and (65), respectively, given below:

{tilde over (y)}′=y′(t,y,{tilde over (t)})=y+[v _(T)(t,y)n _(+c)+u(t,y)]({tilde over (t)}−t)+½g({tilde over (t)}−t)²,  (64)

and

{tilde over (v)} ₊′({tilde over (t)})=v ₊′(t,y,{tilde over (t)})=v_(T)(t,y)n _(+c) +u(t,y)+g({tilde over (t)}−t),  (65)

where n_(+c)=±1 is a unit vector of arbitrary direction with the initialpoint of origin y at time t.

The step of determining a traveling time needed, for each particlediverging from non-moving point y at the given time, t, to cross acontrol surface enclosing point y in point y′ of control surfaceincludes a step of resolving the equation of projectile motion, Equation(66) given below, with respect to the ballistic traveling time, φ₊:

y′=y+[v _(T)(t,y)n _(+c) +u(t,y)]φ₊+½g(φ₊)²,  (66)

which is obtained by assigning {tilde over (t)}=t_(a)′ and substitutionof φ₊=t_(a)′−t in Equation (64) and rearrangement of the terms.Traveling time φ₊ is needed to indicate time t_(a)′ at which thedeparting particle has reached point y′ on control surface 1201 or 1202enclosing point y.

The step of defining an appropriate instant unit vector directingthermal velocity component along the diverging ballistic trajectory ofeach particle diverging from the general non-moving point y at time t insome embodiments includes the steps of:

presenting Equation (66) in the following form:

φ₊ v _(T) n _(+c) =y′−y _(a) ^(c),  (67)

where n_(+c) is a unit vector and y_(a) ^(c) (not shown) is a virtualposition of a particle having zero magnitude of an arbitrary or thermalvelocity at point y at time t, which is observed at time t_(a)′, and, informula form, is expressed as:

y _(a) ^(c) =y+u(t,y)φ₊+½g(φ₊)²;  (68)

and deriving n_(+c) from Equation (67), which is given as:

$\begin{matrix}{{n_{+ c} = \frac{y^{\prime} - y_{a}^{c}}{{y^{\prime} - y_{a}^{c}}}},} & (69)\end{matrix}$

where referring to FIG. 12, n_(+c) (not shown) is interpreted as a unitvector of origin y_(a) ^(c) at time t_(a)′, which directs the vector ofthe thermal velocity component and t_(a)′ is the time of positioning inpoint y′, which is on the control surface 1201 or 1202.

The step of defining the velocity vector, v₊, of each particle at themoment of crossing the control surface in point y′ includes representingv₊ by applying Equation (70) given below:

v ₊(t,y,t _(a) ′,y′)=v _(T)(t,y)n _(+c) +u(t,y)+gφ ₊,  (70)

which is obtained by assigning {tilde over (t)}=t_(a)′ and substitutionof t_(a)′−t=φ₊ in Equation (65).

1.2.3.2 Analytical Tool for Computing the Probability of BallisticTraveling Along the Ballistic Trajectory of the Particle Diverging fromthe General Non-Moving at the Given Time

The step of formulating the probability, Q₊, of ballistic travelingalong the trajectory from the point of the divergence y to point y′ onthe control surface 1201 or 1202 includes a step of representingQ₊(t_(a)′,t) by applying Equation (71) given below:

Q ₊(t _(a) ′,t)=Q ₊(0,φ_(i+))=exp(−∫_(t) ^(t) ^(a) ^(′) P _(c)({tildeover (y)}′({tilde over (t)}))v _(rel)({tilde over (y)}′({tilde over(t)}))d{tilde over (t)}),  (71)

where t_(a)′=t+φ₊ is a time of the particle positioning in point y′within the volume of the model gas system, P_(c)({tilde over(y)}′({tilde over (t)})) is the number of particles within a collisiontube placed in the model gas having particle density of n({tilde over(y)}′({tilde over (t)})) at time {tilde over (t)}, v_(rel)({tilde over(y)}′({tilde over (t)})) is an average magnitude of the relativevelocity between the particle traveling along the diverging ballistictrajectory and another nearby particle being within the collision tubealong the diverging ballistic trajectory at time {tilde over (t)}, withthe method of its determining described above, {tilde over (y)}′(t) isdefined by Equation (64).

1.2.3.3 Analytical Tool for Computing the Net Rate of Property EffluxPer Unit Volume from the General Non-Moving at the Given Time

The step of defining the vector field of the property flux, J_(y) ^(Ψ),along each of the ballistic trajectories in point {right arrow over(r)}′ of the trajectory includes a step of representing J_(y) ^(Ψ) byapplying Equation (72) given below:

J _(y) ^(Ψ)+½n(t,y)Q ₊(t _(a) ′,t)v ₊(t,y,t _(a) ′,y′)Ψ_(λ)(t,y,t _(a)′,y′),   (72)

where

$\frac{\partial}{\partial y}$

-   -    is the divergence operator in one-dimensional configuration,    -   t is the given time,    -   y is the position of the starting point of a plurality of        ballistic trajectories of diverging ballistic particles,    -   y′ is the position of the ending point of one of the plurality        of ballistic trajectories of diverging ballistic particles,    -   t_(a)′ is a time of positioning in the ending point,    -   v₊(t,y,t_(a)′,y is a velocity vector of each of the plurality of        diverging ballistic particles at the time of positioning the        ending point, which is defined by Equation (70),    -   Q₊(t_(a)′,t) is the probability of traveling along the ballistic        trajectory of one of the plurality of diverging ballistic        particles, which is defined by Equation (71),    -   Ψ_(λ)(t,y,t_(a)′,y is property value carried by one of the        plurality of diverging ballistic particles at the time of        positioning the ending point, and    -   n(t,y) is particle density in point in the starting point at the        given time.

The step of defining the net rate of property efflux per unit volume,B_(out) ^(Ψ_FS), from a point source positioned in the generalnon-moving point y at the given time t includes a step of representingB_(out) ^(Ψ_FS) by applying Equation (73) given below:

$\begin{matrix}{{{B_{out}^{{\Psi\_}{FS}}\left( {t,y} \right)} = {\left( {\frac{\partial}{\partial y}J_{y}^{\Psi}} \right)_{y^{\prime}\rightarrow y} = {\frac{1}{2}\left\{ {\frac{\partial}{\partial y}\left\lbrack {{n\left( {t,y} \right)}{Q_{+}\left( {t_{a}^{\prime},t} \right)}{v_{+}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}{\Psi_{\lambda}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}} \right\rbrack} \right\}_{y^{\prime}\rightarrow y}}}},} & (73)\end{matrix}$

which is obtained by executing the operator

$\frac{\partial}{\partial y}$

to the vector field of Equation (72) followed by executing limit y′→y,which also leads to the limit y_(a) ^(c)→y.

1.2.4. Defining the Temporal Rate of the Property Change Per Unit Volumein the General Non-Moving at the Given Time

The step S105 of defining the temporal rate of the property change perunit volume, {dot over (A)}^(fl), in the same general non-moving at thesame time includes a step (not shown) of calculating by applyingEquation (74) given below:

$\begin{matrix}{{{\overset{.}{A}}^{fl} = {\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi_{\lambda}\left( {t,y} \right)}} \right\rbrack}},} & (74)\end{matrix}$

where Ψ_(λ) is scalar or vector property possessed by each particle inpoint y at time t; n is particle density in point y at time t.

1.2.5. Analytical Tool for Defining an Integral Property BalanceEquation in the General Non-Moving Point at the Given Time

The step S106 includes establishing an integral property balanceequation for each of one or more properties being transported by theplurality of particles in each of the plurality of non-moving points ina way that, in the general point at the time, a term of the net rate ofproperty influx per unit volume is equated to sum of the first term oftemporal rate of a property change per unit volume and a second term ofthe net rate of property efflux per unit volume.

What is needed to form property balance in the model gas system is torecognize that each point of collisions is both a collector of propertydelivered by converging ballistic particles from the entire model gassystem and a disperser of property into surrounding taken away bydiverging ballistic particles.

FIG. 13 is a block diagram of a word equation of a general propertybalance in a model gas system according to the principles of the method.According to the method, it is constructed in such a way that the netrate of property influx per unit volume from the surrounding model gassystem, B_(in) ^(Ψ_FS), is equated to the temporal rate of propertychange per unit volume in the same point

${\frac{\partial}{\partial t}\left\lbrack {n\; \Psi_{\lambda}} \right\rbrack},$

and the net rate of property efflux per unit volume from the same pointwhich is a source point for the diverging ballistic particles originatedfrom the collisions in this point, each taking away property intrinsicto the model gas near the same point simultaneously and dispersing itinto surrounding model gas, B_(out) ^(Ψ_FS), which is expressed informula form as:

$\begin{matrix}{{B_{i\; n}^{{\Psi\_}{FS}}\left( {y,t} \right)} = {{B_{out}^{{\Psi\_}{FS}}\left( {y,t} \right)} + {{\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi_{\lambda}\left( {t,y} \right)}} \right\rbrack}.}}} & (75)\end{matrix}$

In a steady-state model gas system, the equation above is reduced toEquation (76) given below,

B _(in) ^(Ψ_FS)(y,t)=B _(out) ^(Ψ_FS)(y,t).  (76)

Concerning FIG. 1 and FIG. 1A, steps S103A, S104, and S105 are needed todefine the property balance equation for a model gas system with nogas-solid interfaces. In this, the step S103 of defining the net rate ofproperty influx per unit volume from the surrounding model gas system,B_(in) ^(Ψ_FS), includes a step of representing B_(in) ^(Ψ_FS) byapplying Equation (77) given below:

B _(in) ^(Ψ_FS) =B _(in) ^(Ψ_F)  (77)

where B_(in) ^(Ψ_F) is computed according to step S103A.

In this embodiment, formulation of the integral form of the propertybalance equation in a one-dimensional model gas system with noboundaries includes applying Equation (78) given below, which isobtained by substitution of Equations (38) and (73) in Equation (75):

$\begin{matrix}\begin{matrix}{{{{- \frac{1}{2}}\frac{\partial}{\partial y}{\int_{y_{b\; 2}}^{y_{b\; 1}}{{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}{\Psi_{i\; n}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{dy}^{\prime}}}} = {{\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi_{\lambda}\left( {t,y} \right)}} \right\rbrack} + {\frac{1}{2}\left\{ {\frac{\partial}{\partial y}\left\lbrack {{n\left( {t,y} \right)}{Q_{+}\left( {t_{a}^{\prime},t} \right)}{v_{+}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}{\Psi_{\lambda}\left( {t,y} \right)}} \right\rbrack} \right\}_{y^{\prime} = y}}}},} & (78)\end{matrix} & \;\end{matrix}$

where in the left-hand of the equation, Q_(i)(t,t_(i)′),Z_(V)(t_(i)′,y′), v(t_(i)′,y′,t,y) are obtained from Equations (21),(32), and (16), respectively, and, in the second right-hand term,Q₊(t_(a)′,t) and v₊(t,y,t_(a)′,y′) are expressed by Equations (71) and(70), respectively.

In the first term of the equation above, the probability of ballistictraveling a distance from point y′ to point y is a limiting factordefining the range of the property influx impact in the propertybalance. When executing numerical calculations by any suitable computersystem, the values of the probability of ballistic traveling are storedand executed if they fall above the minimum absolute value (MAV)corresponding to a suitable format such as Single, Double, or Extended.Symbolically, the above is represented:

Q _(i)(t,t _(i)′)=exp(−∫_(t) _(i) _(′) ^(t) P _(c)(y({tilde over (t)},t_(i)′))v _(rel)(y({tilde over (t)},t _(i)′))d{tilde over(t)})≥MAV.  (79)

The corresponding estimate of the non-equality above is given by:

$\begin{matrix}{{\phi_{i}^{\max} \leq \frac{- {\ln ({MAV})}}{{P_{c}\left( {y,t} \right)}{v_{rel}\left( {y,t} \right)}}},} & (80)\end{matrix}$

where φ_(i_max) is an approximate value of the longest traveling timefor the particle to reach point y at time t and being detected in thenumerical calculation by the computer system.

Also, we have recognized that a predetermined dynamic history duringsome period immediately preceding the initial time t₀ (“past”) willmaximize the accuracy of prediction of the dynamic evolution of themodel gas system especially in the most recent times following theinitial time t₀ (“future”). Under the method, an approximate reliableperiod of time, φ_(i0) ^(rel), with an established dynamic history,which has immediately preceded the initial time t₀, is calculatedaccording to Equation (81):

$\begin{matrix}{{\phi_{i\; 0}^{rel} = {- \frac{\ln ({MAV})}{\left\lbrack {P_{c}v_{rel}} \right\rbrack_{\max}}}},} & (81)\end{matrix}$

where φ_(i0) ^(rel) is the reliable period before the initial time, MAVis a minimum absolute value corresponding to a precision formatsupported by computer hardware, P_(c) is a number of particles within acollision tube of a unit length placed in the model gas, v_(rel) is anaverage magnitude of the velocity of a traveling particle with respectto one of the plurality of nearby particles, and [P_(c)v_(rel)]_(max) isthe highest value of the product of P_(c)v_(rel) in the model gassystem.

Here we have recognized that if a period with a pre-established dynamichistory is less than the approximate reliable period then not allpossible property flux or information from the “past” (t<t₀) will bereliably transported in the “future” (t≥t₀). Therefore, the period ofthe pre-established dynamic history is set to be not less than theapproximate reliable period of time, φ_(i0) ^(rel).

Adding to steps S103A, S104, and S105 of steps S103B and S103D willprovide defining the property balance equation for the model gas systemwith gas-solid interfaces revealing a purely diffuse scattering of theparticles in the “open channel” geometry.

FIG. 6 shows schematically a one-dimensional model gas system, in whichthe main control volume of the unit length volume dV=ΔxΔy is at y withinthe model gas flow confined by two parallel plates at y_(b1) and y_(b2)and having normal vectors n_(b1) and n_(b2), respectively, that aredirected toward the confined model gas. Here we assume diffusivescattering with no adsorption effects after colliding with the surface.

For further clarity, the property balance of model gas flow inone-dimensional configuration with a lack of the external field of forcein the y-direction (B_(in) ^(Ψg_FS)=0) is analyzed.

In this, the step S103 of defining the net rate of property influx perunit volume from the surrounding model gas system, B_(in) ^(Ψ_FS),includes, for example, a step of representing B_(in) ^(Ψ_FS) by applyingEquation (82) given below:

B _(in) ^(Ψ_FS) =B _(in) ^(Ψλ_F) +B _(in) ^(Ψλ_bd1) +B _(in) ^(Ψλ_bd2)+{dot over (Q)} _(Ψ),  (82)

where B_(in) ^(Ψλ_F), B_(in) ^(Ψλ_bd1)+B_(in) ^(Ψλ_bd2), and {dot over(Q)}_(Ψ) are the net rate of property influxes calculated by performingsteps S103A, S103B, and S103D, respectively.

In this embodiment, formulation of the integral form of the propertybalance equation in a one-dimensional model gas system having boundarieswith diffuse scattering from model gas-surface interfaces includesapplying Equation (83) given below, which is obtained from Equation (75)by substitution of it's the first left-hand term by Equation (82):

$\begin{matrix}{{{{B_{in}^{\Psi \; {{\lambda\_}F}}\left( {y,t} \right)} + {B_{in}^{\Psi \; {{\lambda\_}b1}}\left( {y,t} \right)} + {B_{in}^{\Psi \; {{\lambda\_}b2}}\left( {y,t} \right)} + {{\overset{.}{Q}}_{\Psi}\left( {y,t} \right)}} = {{B_{out}^{\Psi \; {{\lambda\_}{FS}}}\left( {y,t} \right)} + {\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi_{\lambda}\left( {t,y} \right)}} \right\rbrack}}},} & (83)\end{matrix}$

which further expands by substitution of the first, second, and thethird left-hand terms by Equations (38), (49), and (50), respectively,and by substitution of the first right-hand terms by Equation (73):

$\begin{matrix}{{{{- \frac{1}{2}}\frac{\partial}{\partial y}{\int_{y_{b\; 2}}^{y_{b\; 1}}{{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}{\Psi_{in}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}\ {dy}^{\prime}}}} - {\frac{\partial}{\partial y}\left\lbrack {{Q_{ib}\left( {t,t_{{ib}\; 1}^{\prime}} \right)}{Z_{b}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}\frac{v_{b\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1},t,y} \right)}{v_{Tb}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1}} \right)}{\Psi_{\lambda \; b\; 1}\left( {t_{{ib}\; 1}^{\prime},y_{b\; 1},\lambda,\phi_{{ib}\; 1}} \right)}} \right\rbrack} - {\frac{\partial}{\partial y}\left\lbrack {{Q_{ib}\left( {t,t_{{ib}\; 2}^{\prime}} \right)}{Z_{b}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2}} \right)}\frac{v_{b\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2},t,y} \right)}{v_{Tb}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2}} \right)}{\Psi_{\lambda \; b\; 2}\left( {t_{{ib}\; 2}^{\prime},y_{b\; 2},\lambda,\phi_{{ib}\; 2}} \right)}} \right\rbrack} + {{\overset{.}{Q}}_{\Psi}\left( {y,t} \right)}} = {{\frac{1}{2}\left\{ {\frac{\partial}{\partial y}\left\lbrack {{n\left( {t,y} \right)}{Q_{+}\left( {t_{a}^{\prime},t} \right)}{v_{+}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}{\Psi_{\lambda}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}} \right\rbrack} \right\}_{y^{\prime}\rightarrow y}} + {\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi \left( {t,y} \right)}} \right\rbrack}}} & (84)\end{matrix}$

1.3 the Fluid Model of Fluids in a Stable Condition

We have also recognized that the model gas system can be in a stablecondition or a steady state, in which properties or variablescharacterizing the model gas system and defining its behavior areunchanged in time. It implies that properties or variablescharacterizing the model gas system are dependent on the position inspace but not time. We have recognized the model gas system may be in asteady state before the initial time and become a temporal model gassystem at times after the initial time if it is affected by an externalforce or geometry modification at the moment of the initial time.Approximating fluid by the model gas in a stable condition in someembodiments further includes:

treating each point of the plurality of points in space occupied by themodel gas as a point of collisions for each of the plurality ofparticles having the ballistic trajectory having the ending point in oneof the plurality of points;

treating each of a plurality of the points of collisions as either apoint source for each of a plurality of diverging particles or a pointsink for each of a plurality of converging particles; and

treating each particle of the plurality of particles moving from thepoint source to the point sink as a property carrier created in thepoint source by obtaining one or more of properties of specific valuesbeing intrinsic to the model gas surrounding the point source, and endedin the point sink by transferring one or more of properties of specificvalues in the point sink,

where each of the plurality of particles is assigned to travel with aprobability between any of two points of a plurality of points ofconsecutive collisions in space occupied by the model gas by following aballistic trajectory governed by a law of motion in free space, theballistic trajectory having a starting point in one of a plurality ofpoints of original collisions and an ending point in one of a pluralityof points of ending collisions, and

where the velocity of each of a plurality of point sources is assignedto be equal to the mass flow velocity of the model gas in each of aplurality of corresponding points of original collisions;

treating each of a plurality of collisions on a gas-solid interface ofthe model gas system, which has resulted in diffuse particle scatteringfrom the gas-solid interface, as an act of interaction involving aproperty transfer from the gas-solid interface to a scattered particle;and

treating each of a plurality of points of diffuse particle scattering onthe gas-solid interface as a heterogeneous point source for each of aplurality of scattered particles,

where the gas-solid interface reveals mixed diffuse and specularparticle scatterings,

where the velocity of each of a plurality of heterogeneous point sourceson the gas-solid interface is assigned to be equal to the velocity ofthe gas-solid interface in each of a plurality of corresponding pointsof diffuse particle scattering, and

where point source strength of each of the plurality of heterogeneouspoint sources on the gas-solid interface is assigned to be directlyproportional to a property accommodation coefficient in each of theplurality of corresponding points of diffuse particle scattering, theproperty accommodation coefficient which is probability, for an incidentparticle, to accommodate one or more of properties intrinsic to thegas-solid interface and to scatter back in the model gas as a diffuseparticle, the property accommodation coefficient being in a range fromzero to one.

1.4 the Model of the Property Balance in a Stable Condition

In some embodiments, modeling transport processes for model gas flow ina stable condition includes the steps of:

specifying geometry model and boundary conditions;

defining a net rate of property influx per unit volume from the modelgas system in a general point of a plurality of non-moving points inspace occupied by the model gas, the net rate of property influx perunit volume, which is formed by way of a plurality of convergingballistic particles, each of the plurality of converging ballisticparticles is selected from the plurality of particles by the ballistictrajectory having the ending point in the general point of the pluralityof non-moving points in space occupied by the model gas,

defining a net rate of property efflux per unit volume from the generalpoint, the net rate of property efflux per unit volume, which is formedby way of a plurality of diverging particles, each of the plurality ofdiverging ballistic particles is selected from the plurality ofparticles by the ballistic trajectory having the starting point in thegeneral point of the plurality of non-moving points;

forming an integral property balance equation for each of one or moreproperties being transported by the plurality of particles in each ofthe plurality of non-moving points,

where the integral property balance equation is established in a waythat, in the general point, a term of the net rate of property influxper unit volume is equated a term of the net rate of property efflux perunit volume; and

computing the flow of the model gas,

where computing the flow of the model gas includes resolving a system ofone or more of integral property balance equations in each of theplurality of non-moving points.

In the stable model gas system not affected by the external field offorce, with model gas confined between two infinite parallel plateshaving gas-solid interfaces with identical mixed diffuse and specularscattering components, formulating the property balance includesapplying Equation (85) given below:

$\begin{matrix}{{{B_{in}^{{\Psi\_}F} + B_{in}^{{{\Psi 0\_}F}{\_ {bs}1}} + B_{in}^{{{\Psi 0\_}F}{\_ {bs}2}} + {{\overset{.}{Q}}_{\Psi}\left( {t,y} \right)}} = {\frac{1}{2}\left\{ {\frac{\partial}{\partial y}\left\lbrack {{{n\left( {t,y} \right)}Q} + {\left( {y,y^{\prime}} \right)v} + {\left( {y,y^{\prime}} \right){\Psi \left( {y,y^{\prime}} \right)}}} \right\rbrack} \right\}_{y^{\prime}\rightarrow y}}},} & (85)\end{matrix}$

which is obtained by substitution of Equations (73) and (59) in Equation(75) and considering that

$\begin{matrix}{{\frac{\partial}{\partial t}\left\lbrack {{n\left( {t,y} \right)}{\Psi \left( {t,y} \right)}} \right\rbrack} = 0} & (86)\end{matrix}$

The specific embodiments of the method for modeling transport processesin fluids are disclosed more fully in the following model of thecomputer for simulating the fluid flow and non-limiting examples(Simulation 1 and Simulation 2) of the best mode of carrying out theinvention.

1.5 the Model of the Computer for Simulating the Fluid Flow

The present analytical tools and method are further described by way ofthe following one-dimensional fluid flow computer simulations, which areprovided by way of illustration and are not intended to limit thepresent invention.

The computer for simulating the fluid flow includes a computer readablemedium for storing data associated with a computer program, and aprocessor for executing a sequence of instructions from the computerprogram, in which the computer program is configured to use the fluidmodel and the model of property balance.

These embodiments are limited to a one-dimensional model gas system,including incompressible model gas at a uniform temperature with noexternal field of the force applied in y direction. The presentinvention, however, is not limited to the case, but can be preferablyapplied to two- and three-dimensional configurations of the model gassystem with no limitations mentioned above. The method can also begeneralized to the model gas flow affected by the external field offorce with no difficulty.

In these simulations and computations, consideration has been made thattransported properties that are not vulnerable to any aging process. Byother words, the property acquired by each particle in a point duringthe preceding collision is delivered unchangeable in the point of thefollowing collision. However, the method can be generalized to thetransport of property/properties having limited life expectancy with nodifficulty.

FIG. 6 illustrates the diversity of ballistic trajectories of particlesinvolving in the property transport by the present method in theone-dimensional configuration. Ballistic trajectories 603 and 606 showmovement of particles in a positive direction. Ballistic trajectories601 and 605 show movement of the particles in the negative direction.

Ballistic trajectories 601 and 603 show movement of particles from themodel gas volume into targeting point, y. Ballistic trajectories 605 and606 show movement of particles from plate 615 and plate 616,respectively, into targeting point y.

Usually, the magnitude of the thermal velocity of the particles is muchhigher than the magnitude of mass flow velocity vector along they-direction, |u|<<v_(T).

The following Table 2, referring to FIG. 6, summarizes symbols and theirmeanings, which are applicable to the simulations and calculationsdescribed below with additional reference to the figures and thespecification.

TABLE 2 Symbol Value Short description v_(T) constant uniformtemperature; thermal velocity is uniform within the model gas system nconstant particle density σ_(c) constant the cross-section of collisionsλ 0 aging or decay coefficient u 0 mass flow velocity along y axis n_(i)$\frac{y - y^{\prime}}{{y - y^{\prime}}}$ arbitrary unit vector withthe point of origin y′ v v_(T)n_(i) the velocity of the particle inpoint y′ n₊ $\frac{y^{\prime} - y}{{y^{\prime} - y}}$ arbitrary unitvector with the point of origin y v₊ the velocity of the particlereaching point y′ after its original collision in point y at time tv_(rel) v_(T) the average magnitude of the relative velocity at theuniform temperature Z_(V) = Z_(V′) ½nP_(c)v_(rel) the rate of collisionsper unit volume per unit time or the particle flux production strengthP_(c) σ_(c)n the number of particles within a collision tube of a unitlength Q_(i) exp(−P_(c)|y−y′|) the probability of ballistic travelingbetween y′ and y Q₊ exp(−P_(c)|y′−y|) the probability of ballistictraveling between y and y′ y_(b1) = H the position vector of plate 615y_(b2) = 0   the position vector of plate 616 n_(b1) = −1 the normalvector to plate 615 n_(b2) = +1 the normal vector to plate 616 n_(ib1) =−1 the direction vector of the particle scattered from plate 615 n_(ib2)= +1 the direction vector of the particle scattered from plate 616v_(Tb1), v_(Tb2) v_(Tb1) = v_(Tb2) = v_(T) the magnitude of the velocityof particles scattered from being at rest plate 615 and plate 616 at theuniform temperature v_(b1) v_(T)n_(ib1) the velocity vector of theparticle scattered from plate 615 v_(b2) v_(T)n_(ib2) the velocityvector of the particle scattered from plate 616 Z_(b1) = Z_(b2)$\frac{Z_{V}}{2\; P_{c}}$ the rate of collisions on plate 615 or plate616 at the uniform temperature Q_(ib1) exp(−P_(c)(H − y)) theprobability of ballistic traveling from y_(b1) = H to y Q_(ib2)exp(−P_(c)y) the probability of ballistic traveling from y_(b2) = 0 to yt′_(ib2) $t - \frac{y}{v_{T}}$ the time of scattering from plate 616 forthe particle targeting point y and time t t′_(i)$t - \frac{{y - y^{\prime}}}{v_{T}}$ the time of an original collisionin point y φ′_(i) $\frac{{y - y^{\prime}}}{v_{T}}$ a traveling periodneeded to travel from y′ to y φ′_(ar) $\frac{y^{\prime}}{v_{T}}$ atraveling period needed for the converging ballistic particle to deliverproperty from plate 616 to point y′ φ_(ar) $\frac{y}{v_{T}}$ a travelingperiod needed for the converging ballistic particle to deliver propertyfrom plate 616 to point y Ψ_(0b2) the value of property obtained byadiffusely scattered particle from plate 616.

To formulate a general integral form of the property balance adapted toa particular model gas system in a general non-moving point y at thegiven time, t, there is need to perform the following.

First, using 1.2.1.1 “Analytical tool for computing the trajectory andvelocity of a ballistic movement . . . ,” a trajectory and associatedwith it both an appropriate instant unit vector directing thermalvelocity component, the point of the original collision in y′, the time,t_(i)′, of the original collision are formulated. This is done by thesteps of:

(1) formulating the velocity of the particle along y-direction arrivingin point y at time t is calculated from Equations (87) given below,which is obtained from Equation (16):

v(t,y)=v _(T) n _(i),  (87)

where n_(i) is a unit vector of arbitrary direction and is expressed:

$\begin{matrix}{{n_{i} = \frac{y - y^{\prime}}{{y - y^{\prime}}}},} & (88)\end{matrix}$

(2) determining the position of the particle at traveling time {tildeover (t)} is calculated from Equation (89) given below, which isobtained from Equation (10):

{tilde over (y)}=y′+v _(T)({tilde over (t)}−t _(i)′)n _(i) =v _(T){tilde over (t)}+y′−v _(T) t _(i)′  (89)

for the particles traveling in the positive direction and

{tilde over (y)}=y′−v _(T)({tilde over (t)}−t _(i)′)=−v _(T) {tilde over(t)}+y′+v _(T) t _(i)′  (90)

for the particles traveling in the negative direction;

(3) determining ballistic traveling time from point y′ to point y,φ_(i), by applying Equation (91) given below, which is obtained fromEquation (12) in which the impacts of g and u are neglected:

$\begin{matrix}{{\phi_{i} = \frac{{y - y^{\prime}}}{v_{T}}};} & (91)\end{matrix}$

(4) finding the departing time of the particle at y′, t_(i)′, which iscalculated from Equation (92) given below:

t _(i) ′=t−φ _(i).  (92)

Second, using 1.2.1.2 “Analytical tool for computing the probability ofballistic traveling . . . ,” probability function quantifying theprobability of the ballistic traveling along the ballistic trajectory isformulated.

Here, upon recognizing the magnitude of the relative velocity betweenthe component of the mass flow velocity of the traveling and one ofnearby particles is insignificant in comparison with the magnitudes ofthermal velocities of the particles, we obtain:

(1) the average magnitude of relative velocity between the particles atthe uniform temperature is calculated from Equation (93) given below,which is obtained from Equation (30):

v _(rel) =v _(T);  (93)

(2) the survival probability that a particle will have traveled alongthe ballistic trajectory in incompressible model gas at the uniformtemperature from y′ and y is calculated from Equation (94) given below,which is obtained by substitution of Equation (93) in Equation (21) andconsidering that P_(c)=constant:

Q _(i)(y,y′)=exp(−P _(c) |y−y′|);  (94)

(3) the survival probability, Q_(ib1)(y,H), that a particle will havetraveled along the ballistic trajectory in incompressible model gas atthe uniform temperature from y_(b1)=H and y is calculated from Equation(95) given below:

Q _(ib)(t,t _(ib1)′)=Q _(ib1)(y,H)=exp(−P _(c)(H−y));  (95)

(4) the survival probability, Q_(ib2)(y,0), that a particle will havetraveled along the ballistic trajectory in incompressible model gas atthe uniform temperature from y_(b2)=0 and y is calculated from Equation(95) given below:

Q _(ib)(t,t _(ib2)′)=Q _(ib2)(y,0)=exp(−P _(c) y).  (96)

Third, using 1.1.3 The model of the geometry model and the dynamichistory, a functional representation of property being transported byballistic particles is adapted to the particular model gas system, whichis characterized by uniform temperature and lack of external field offorce and internal aging processes. Equation (7) is modified as:

Ψ_(in)(t _(i) ′,y′,t,y)=Ψ₀(t _(i) ′,y′).  (97)

For clarity, the value of the property on the upper plate 615, is chosento be

Ψ_(0b1)=0;  (98)

The value of the property on the lower plate 616 is defined as:

Ψ_(0b2)≡=Ψ_(0b2−)[t _(ib2) ′<t ₀]+Ψ_(0b2+)[t _(ib2) ′≥t ₀],  (99)

where t₀ is initial time for the model gas system, which is associatedwith a time of modification of property Ψ_(0b2) on the gas-solidinterface of plate 616. We have recognized that scattered from plate 616diffuse particles will travel and initiate a modification of propertyΨ₀(t_(i)′,y′) upon arriving in point y′. Symbolically, this is expressedaccording to Equation (100), which is obtained by substitution oft_(i0)=t₀+φ_(ar)′ in Equation (9):

Ψ₀(t _(i) ′,y′)=Ψ⁰⁻(t _(i) ′,y′)[t _(i) ′<t ₀+φ_(ar)′]+Ψ(t _(i) ′,y′)[t_(i) ′≥t ₀+φ_(ar)′].  (100)

Finally, using 1.2.5 “Analytical tool for defining an integral propertybalance equation in the general non-moving point at the given time,” amodified integral form of the property-balance equation, adapted to aparticular model gas system, is formulated. Specifically, forincompressible model gas expanding in y-direction between parallelinfinite plates being spaced at a distance H in the open channelgeometry at the uniform temperature the modified integral for of theproperty-balance equation is formulated by Equation (101) given below,which is obtained by substitution of symbols or corresponding values ofTable 2 and of Equations (87)-(100) in Equation (84), then executingdifferentiation while considering that y′≠y:

$\begin{matrix}{{{n{\frac{\partial}{\partial t}\left\lbrack {\Psi \left( {t,y} \right)} \right\rbrack}} + {Z_{V}{\Psi \left( {t,y} \right)}}} = {{{- \frac{1}{2}}Z_{V}\frac{\partial}{\partial y}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{\Psi_{0 -}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} < {t_{0} + \phi_{ar}^{\prime}}} \right\rbrack}{dy}^{\prime}}}} + {\frac{Z_{V}}{2}{{\exp \left( {{- P_{c}}y} \right)}\left\lbrack {{\Psi_{{0b\; 2} -}\left\lbrack {t_{{ib}\; 2}^{\prime} < t_{0}} \right\rbrack} + {\Psi_{{0b\; 2} +}\left\lbrack {t_{{ib}\; 2}^{\prime} \geq t_{0}} \right\rbrack}} \right\rbrack}} - {\frac{1}{2}Z_{V}\frac{\partial}{\partial y}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{\Psi \left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} \geq {t_{0} + \phi_{ar}^{\prime}}} \right\rbrack}{{dy}^{\prime}.}}}}}} & (101)\end{matrix}$

In the equation above, the first right-hand term quantifies apre-established rate of the property influx formed by convergingballistic particles, each delivering property obtained from acorresponding point source being advanced to the initial time. Thesecond right-hand term quantifies a pre-established rate of the propertyinflux formed by converging ballistic particles, each deliveringproperty obtained from a corresponding heterogeneous point source on thegas-solid interface of plate 616. The third right-hand term formulates arate of the property influx formed by converging ballistic particles,each delivering property obtained from a corresponding point sourcebeing the same as or later than the initial time.

In some embodiments, the step of forming the integral property balanceequation further includes an approximating step of approximating theintegral property balance equation by neglecting the first term of thetemporal rate of the property change.

Specifically, by if, in the general non-moving point at the given time,a temporal rate of the property change per unit volume is negligiblecompared to the rate of collisions in the general non-moving point,which is true in majority of realistic conditions in the model gassystem, then Equation (101) is reduced to:

$\begin{matrix}{{\Psi \left( {t,y} \right)} \cong {{\frac{1}{2}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{\Psi_{0 -}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} < {t_{0} + \phi_{ar}^{\prime}}} \right\rbrack}{dy}^{\prime}}}} + {\frac{1}{2}{{\exp \left( {{- P_{c}}y} \right)}\left\lbrack {{\Psi_{{0b\; 2} -}\left\lbrack {t_{{ib}\; 2}^{\prime} < t_{0}} \right\rbrack} + {\Psi_{{0b\; 2} +}\left\lbrack {t_{{ib}\; 2}^{\prime} \geq t_{0}} \right\rbrack}} \right\rbrack}} + {\frac{1}{2}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{\Psi \left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} \geq {t_{0} + \phi_{ar}^{\prime}}} \right\rbrack}{{dy}^{\prime}.}}}}}} & (102)\end{matrix}$

A step of formulating a solution for the first approximation of amodified integral form of the property-balance equation includes stepsof:

(1) modify Equation (101) by assignment, in the last right-hand term of

Ψ=0,  (103)

which is a zero approximation for a solution of Equation (101);

(2) resolve the modified Equation (101) and obtain a solution Ψ₁(t,y)for each of a plurality of points in space occupied by the model gas,which quantifies established effects of external fields capable ofmodifying the property content, prehistoric conditions of the model gasflow, and an established effect of an available gas-solid interface.Based on Equation (102), the first approximation, Ψ₁(t,y), has a form:

Ψ₁(t,y)≅½P _(c)∫₀ ^(H) exp(−P _(c) |y−y′|)Ψ⁰⁻(t _(i) ′,y′)[t _(i) ′<t₀+φ_(ar)′]dy′+½ exp(−P _(c) y)[Ψ_(0b2−)[t _(ib2) ′<t ₀]+Ψ_(0b2+)[t_(ib2) ′≥t ₀]].  (104)

A step of formulating a solution for a k^(th) approximation of amodified form of the property-balance equation includes presenting anapproximated solution for the k^(th) approximation, Ψ_(k)(t,y), for eachof the plurality of points of the space occupied by the model gas insuch form that the left-hand of an equation defines k^(th) approximationthat Ψ_(k)(t,y) and the right-hand of the equation is defined by(k−1)^(th) approximation, which has been obtained in the preceding cycleof the approximation, which is explained symbolically by Equation (105)given below:

Ψ_(k)(y,t)=Ψ₁(y,t)+Ψ_(ak)(y,t),  (105)

where Ψ_(ak)(y,t) is an adjustment to the first approximation ofvelocity Ψ₁(y,t), which is associated with the mutual effect of theproperty exchange between particles surrounding point y at time t, whichis expressed as follows

Ψ_(ak)(y,t)=½P _(c)∫₀ ^(H) exp(−P _(c) |y−y′|)Ψ_(k−1)(t _(i) ′,y′)[t_(i) ′≥t ₀+φ_(ar)′]dy′.  (106)

The property value in a specific point of space occupied by the modelgas at the given time results from the impact of both the rate ofproperty influx formed by particles being directly scattered fromconfining one or more of gas-solid interfaces and the rate of propertyinflux formed by particles that have acquired the property from each oftheir preceding collision being within the space occupied by the modelgas.

2. Performing Flow Computation by Using the Computer for Simulating theFluid Flow

The step S107 of performing flow calculation is further embodied inthese two numerical simulations, which are provided for the illustrationof the method and are construed not to limit the scopes of theinvention. Simulation 1 and Simulation 2 of the following sectionprovide detailed steps of computation of velocity profiles in theone-dimensional incompressible flow. However, the method can begeneralized to a three-dimensional model gas with no difficulty.

Embodiments within the scope of the present invention also includecomputer program means including the computer readable medium havingrepresented in that program code means. Program code means includesexecutable instructions and data that cause a computer to perform acertain function of a group of functions either directly or byconversion to another language, including reproduction in a differentmaterial form. All these simulations suggest the applicability of theinnovative analytical tools to construct specific Computational FluidDynamics models.

A method for modeling transport processes and computerized simulatingthe flow of the model gas includes, for example, steps of:

establishing an initial time,

where the initial time is treated as a time of specific modification ofthe model gas system in a specific location including specificmodification of one or more of temperature, the velocity of thegas-solid interface, and application of an external field that modifiesproperty content;

specifying the given time,

where the given time is greater than or equal to the initial time;

specifying the model gas system,

where the model gas system is specified by defining model gas propertiesincluding gas material, pressure, temperature and by establishing ageometry model,

where establishing the geometry model includes setting geometry andboundary conditions during a period from a pre-initial time until thegiven time;

establishing discretization parameters including a set of discretizationpoints in space occupied by the model gas,

where the set of discretization points includes a set of non-movingpoints;

calculating, for each of the set of discretization points in spaceoccupied by the model gas, a local initial time,

where the local initial time is greater than or equal to the initialtime;

obtaining, for each of the set of discretization points in spaceoccupied by the model gas at a given advanced time, which is ahead ofthe local initial time, a past property content,

where the past property value is obtained by calculations or bymeasurements;

formulating, for each of the set of discretization points at the giventime, an approximated integral form of property balance,

where the approximated integral form of property balance is establishedin a way that, in a selected point of the set of discretization pointsat the given time, a net rate of property influx per unit volume fromthe model gas system, which is formed by a set of converging particlesand calculated by summing a rate of property influx from the set ofdiscretization points surrounding the selected point is equated to a netrate of property efflux per unit volume from the selected point of theset of discretization points,

where the set of converging particles converging in the selected pointof the set of discretization points at the given time includes a set ofconverging ballistic particles and a second set of converging ballisticparticles, each of the set of converging ballistic particles isdelivered one or more of past property contents from the model gassystem and each of the second set of converging ballistic particles isdelivered one or more of present property contents,

where each of present property contents, Ψ, is obtained from acorresponding point of original collisions at a time, which is greaterthan or equal to the local initial time,

where Ψ is treated as a variable; and

resolving with respect to the variable, for each of the set ofdiscretization points in space occupied by the model gas, theapproximated integral form of property balance and calculating the valueof the transported property,

where resolving and calculating is performed by an approximation method,including a successive approximation method.

Specific examples of the execution of velocity profile in model gasusing Microsoft Visual Basic for Applications with Excel and computed ona PC are shown in FIG. 15 and FIG. 16.

FIG. 14 is an example of a flowchart of the flow calculation by thenumerical method. The numerical method of flow calculation includes theapplication of the computer and is applicable to both Simulation 1 andSimulation 2. The numerical calculation, for example, includes steps of:

specifying the model gas system (S1401);

forming u_(x)-momentum balance equation (S1402);

formulating a solution for the first approximation (S1403);

formulating a solution for the k^(th) approximations (S1404);

establishing discretization parameters (S1405);

providing an error estimation rule (S1406);

computing and storing the values for the first approximation u_(x1s) forcorresponding Ys as an array (S1407):

executing computation by the sequential approximations (S1408); and

displaying the results (S1409).

2.1 Example 1: Simulation 1

In the first simulation, Simulation 1, transient plane Couette flow iscalculated.

FIG. 6 is a schematic view in section of a Newtonian model gas betweenparallel infinite plates in the open channel at uniform temperature andno external field of force applied in y direction and a lack of internalaging processes, namely, λ=0, each particle conserves propertiesacquired by it at a point and time of the original collision and deliverthat properties unchanged to another location upon its final collisionwith another particle in a point sink. The following simulation is shownto illustrate the inventive approach in determining velocity profilesu_(x)(y,t) or velocity distribution across the space between the plates.

A simulation is performed for incompressible model gas expanding in they-direction between parallel infinite plates being spaced at thedistance H in the open channel at the uniform temperature where the topand bottom plates and the model gas are initially at rest. Then, at timet₀=0 the bottom plate starts a sudden movement in x-direction along itsplane with constant speed v_(x0). The velocity profile induced in themodel gas due to the sudden motion of the bottom plate 616 is to bedetermined. Referring to FIG. 14, more detailed description of the stepsof performing flow calculation by the computer is given below.

Step S1401 of specifying the model gas system, for example, is given:

Referring again to FIG. 6, which shows a schematic view in section of aNewtonian model gas between parallel infinite plates in the open channelat a uniform temperature, a case when the top and bottom plates and theliquid are initially at rest is considered. For simplicity, theparticles are considered to have a unit mass.

These parameters describing the model gas system are specified:

kind of model gas (providing a mass of molecules, effectivecross-section of the collisions);

T, pressure P;

the distance between parallel plates, H;

the magnitude and direction of the velocity of the bottom plate inx-direction, which begins the movement at time t≥0 with constantvelocity v_(x0); and allowable error of computation, Max_Er.

The upper plate is at rest, which implies that the velocity of the upperplate in x-direction is zero, namely, v_(xH)=0. At time t<0 theincompressible model gas had uniform temperature T and pressure P.

Step S1402 of forming u_(x)-momentum balance equation includes, forexample, steps of:

(1) defining u_(x0) momentum, which is acquired by each particle in apoint during its preceding collision and transported into a point of thefollowing collision at the given includes applying Equation (107) givenbelow, which is obtained by substitution in Equation (100) of Ψ₀≡u_(x0),Ψ⁰⁻≡u_(x0−)=0, a predetermined initial or prehistoric x-momentumin each point y′ at the time of the departure t_(i)′<0 and Ψ≡u_(x), thex-momentum in each point y′ at the time of the departure t_(i)≡t_(i)′≥0, which is a function to be determined:

u _(x0)(t _(i) ′,y′)=u _(x)(t _(i) ′,y′)[t _(i)′≥0];  (107)

(2) defining momentum, which is acquired by each particle in a pointduring its preceding diffuse scattering from the gas-solid interface ofthe confining plate 615 and transported into a point of the followingcollision at the given time includes applying Equation (108) givenbelow:

Ψ_(0b2) ≡=v _(x0)[t _(ib2)′≥0].  (108)

The upper plate 615 is at rest all the time, which implies that

Ψ_(0b1) ≡v _(xH)≡0;  (109)

(3) forming an integral form of the u_(x) momentum balance equationapplied to a specific point in space at the given time includes applyingEquation (110) given below, which is obtained by substitution inEquation (101) of {dot over (Q)}_(Ψ)=0, and Ψ₀≡u_(x)(t_(i)′,y′)[t_(i)′≥0], Ψ_(0b1)=0, and Ψ_(0b2)=v_(x0)[t_(ib2)′≥0]in the right-hand of the equation and to Ψ₀≡u_(x), in the left-hand ofthe equation:

$\begin{matrix}{{{n\frac{\partial}{\partial t}{u_{x}\left( {t,y} \right)}} + {Z_{V}{u_{x}\left( {t,y} \right)}}} = {{{+ \frac{Z_{V}}{2}}v_{x\; 0}{{\exp \left( {{- P_{c}}y} \right)}\left\lbrack {t_{{ib}\; 2}^{\prime} \geq 0} \right\rbrack}} + {\frac{1}{2}Z_{V}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{u_{x}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} \geq t_{0}} \right\rbrack}{{dy}^{\prime}.}}}}}} & (110)\end{matrix}$

The equation above is further reduced by executing some algebraoperation, which results in:

$\begin{matrix}{{{{\frac{\partial}{\partial t}u_{x}} + {\frac{Z_{V}}{n}u_{x}}} = {{\frac{Z_{V}}{2n}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{u_{x}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {t_{i}^{\prime} \geq 0} \right\rbrack}{dy}^{\prime}}}} + {g_{0}\left( {y,t} \right)}}},} & (111)\end{matrix}$

where a is a frequency factor quantified as

$\begin{matrix}{{a = {\frac{Z_{V}}{n} = {\frac{1}{2}P_{c}v_{T}}}},{and}} & (112) \\{{g_{0}\left( {y,t} \right)} = {\frac{a}{2}v_{x\; 0}{\exp \left( {{- P_{c}}y} \right)}{{H\left( {t - \frac{y}{v_{T}}} \right)}.}}} & (113)\end{matrix}$

Step S1403 of formulating solutions for the first approximationincludes, for example, steps of:

(1) assigning in the first right-hand term of Equation (111)

u _(x)=0,  (114)

which is zero approximation for the velocity of each particlesurrounding point y, and which results in a reduction of Equation (111)to

$\begin{matrix}{{{{\frac{\partial}{\partial t}u_{x\; 1}} + {a\; u_{x\; 1}}} = {{{g_{0}\left( {y,t} \right)}\mspace{14mu} {with}\mspace{14mu} {u_{x\; 1}\left( {y,{t = 0}} \right)}} = 0}};} & (115)\end{matrix}$

(2) formulating a solution to the equation above by applying Equation(116):

u _(x1)(y,t)=exp(−at)∫₀ ^(t) exp(as)g ₀(y,s)ds  (116)

which upon integration results in:

$\begin{matrix}{{{u_{x\; 1}\left( {y,t} \right)} = {\frac{1}{2}v_{x\; 0}{{\exp \left( {{- P_{c}}y} \right)}\left\lbrack {{t - \frac{y}{v_{T}}} \geq 0} \right\rbrack}\left\{ {1 - {\exp \left( {{- \frac{1}{2}}P_{c}{v_{T}\left( {t - \frac{y}{v_{T}}} \right)}} \right)}} \right\}}};} & (117)\end{matrix}$

Step S1404 of formulating solutions for the k^(th) approximationincludes, for example, steps of:

(1) modifying the u_(x)-momentum balance equation, which is shown byEquation (111) in such form that u_(x) in the left-hand of the equationdefines k^(th) approximation, u_(xk), and u_(x)(t_(i)′,y′) in theright-hand of the equation defines (k−1)^(th) approximation, u_(x(k−1)),so Equation (111) is modified to

$\begin{matrix}{\mspace{79mu} {{{{\frac{\partial}{\partial t}u_{x\; k}} + {a\; u_{x\; k}}} = {{g_{{({k - 1})} +}\left( {y,t} \right)} + {g_{0}\left( {y,t} \right)}}},\mspace{79mu} {where}}} & (118) \\{{{g_{{({k - 1})} +}\left( {y,t} \right)} = {\frac{a}{2}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{u_{x{({k - 1})}}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {{t_{i}^{\prime} - \frac{y^{\prime}}{v_{T}}} \geq 0} \right\rbrack}\ {dy}^{\prime}}}}},} & (119)\end{matrix}$

where g₀(y,t) is a well-defined function at any y and t, which definesestablished effects of both external fields and prehistoric conditionsof the model gas and model gas flow and g_((k−1)+)(y,t) is a functioncharacterizing the effect from the events of the original collisionstaken place following the initial time;

(2) presenting an approximated solution for the k^(th) approximation,u_(xk), for each point of the space occupied by the model gas byEquation (120) given below:

u _(xk)(y,t)=u _(x1)(y,t)+u _(ak)(y,t),  (120)

where u_(ak)(y,t) is the adjustment to the first approximation ofvelocity u_(x1)(y,t), which is associated with the mutual effect of themomentum exchange between particles surrounding point y at time t, whichis expressed as follows

$\begin{matrix}{{u_{ak}\left( {y,t} \right)} = {\frac{1}{2}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{{u_{{x{({k - 1})}} +}\left( {t_{i}^{\prime},y^{\prime}} \right)}\left\lbrack {{t_{i}^{\prime} - \frac{y^{\prime}}{v_{T}}} \geq 0} \right\rbrack}{{dy}^{\prime}.}}}}} & (121)\end{matrix}$

Equation (121) is obtained by analytical integration. However, numericalintegration of a general function g₀(y,t) over the time can be performedwith no difficulty, as apprehended by those skilled in the art.

Step S1405 of establishing discretization parameters includes, forexample, steps of:

specifying discretization parameters including the dimension of thesystem for calculation [0,H], the number of space intervals S, testingtime t, and the number of time intervals j; and

dividing the space intervals [0,H] into S equal parts of length l, sothat, for y_(s) ≡y′=sl>y, the range is divided as y=y_(s)<y_(s+1)< . . .<y_(S)=H and, for y_(s) ≡y′=sl<y, the range is divided as 0=y₀<y₁< . . .<y_(s)=y. Analogously, the time interval [0,t] is divided into j equalparts of length h: 0=t₀<t₁ . . . <t_(s)< . . . <t_(j)=t. In the above, sis the index of a discretization point positioned in y′.

Step S1406 of providing an error estimation rule includes, for example,a step of providing the error estimation rule for k^(th) cycle of thesequential approximation from Equation (122) given below

$\begin{matrix}{{Er} = \sqrt{\frac{1}{s}{\sum_{1}^{s}\left( \frac{{u_{ak}\left( {y_{s},t} \right)} - {u_{a{({k -})}}\left( {y_{s},t} \right)}}{u_{xk}\left( {y_{s},t} \right)} \right)^{2}}}} & (122)\end{matrix}$

Step S1407 of computing and storing the values u_(x1s)(y_(i),t) forcorresponding y_(s) as an array includes, for example, presentingEquation (121) in a discrete form by computing the values u_(x1s) forcorresponding y_(i) surrounding point:

$\begin{matrix}{u_{x\; 1s} = {\frac{1}{2}v_{x\; 0}{{{{\exp \left( {{- P_{c}}y_{s}} \right)}\left\lbrack {{t_{i}^{\prime} - \frac{y_{s}}{v_{T}}} \geq 0} \right\rbrack}\left\lbrack {1 - {\exp \left( {{- \frac{1}{2}}P_{c}{v_{T}\left( {t_{i}^{\prime} - \frac{y_{s}}{v_{T}}} \right)}} \right)}} \right\rbrack}.}}} & (123)\end{matrix}$

Step 1408 of executing computation for each cycle k of the sequentialapproximations includes, for example, steps of:

assigning u_(xks+)=u_(x1s) (S1408A);

storing the values of u_(xks+) for corresponding y_(s) (S1408B);

computing u_(aks)−u_(ak)(y_(s),t) for corresponding y_(s) according toEquation (121) (S1408C);

computing u_(xks)=u_(aks)+u_(x1s) according to Equation (120) andstoring the values of u_(xks) for corresponding y_(s) (S1408D);

computing Er according to Equation (122) given above (S1408E);

defining a condition for ending approximations (S1408F),

where If Er>Max_Er, assigning u_(xks+)=u_(xks), and go to (S1408B) and

if ≤Max_Er, the cycle of the sequential approximations is stopped, thencopying u_(xks) for corresponding y_(s), and

displaying the results of computation (S1409).

The sign + in the subscript indicates that the function is used for thefollowing integration with y′ as an argument.

FIG. 15 is a chart showing the velocity profiles across between theplates at different moments of time, which are obtained by a computersimulation performed under the steps listed above, and a comparativeanalytical velocity profile, which is obtained independently by usingthe analytical tools disclosed in the specification for simplifiedstable conditions involving incompressible gas at the uniformtemperature. As shown, the result of 160 sequential approximationscomputed according to Simulation 1 for the testing time t=100 mscoincides with the analytical velocity profile. The computer simulationwas made for the nitrogen gas flow in a channel confined by parallelplates spaced at the distance of 1 meter at the uniform temperature of300K and pressure of 0.6 Pa. The computer simulation by the methoddemonstrates the feasibility of the method and shows significantimprovement from the prior art, which consists, first, in the capabilityto model and compute flows in rarefied and dilute gases, second, in theimprovement of time resolution of computations for unsteady flows.Specifically, the results of the numerical calculations demonstrate thecapability of the method to compute the velocity flow profiles at thetesting time as short as one millisecond.

2.2 Example 2: Simulation 2

In the second simulation, Simulation 2, the steady velocity profileacross between the parallel plates with mixed diffuse and specularparticles scatterings from the plates being at rest is calculated.

The particular embodiment is limited to a simulation of incompressiblestable model gas flow confined in the y-direction between parallelinfinite plates at a uniform temperature with no external field of theforce applied in y-direction where the top and bottom plates are atrest. However, the method can be generalized to the model gas flowaffected by the external field of force or kept at non-uniformtemperature with no difficulty, as apprehended by those skilled in theart. The pressure gradient is constant along with the model gas flow inhe x-direction. The velocity profile induced in the model gas due to thepressure gradient is to be determined. Referring to FIG. 14, moredetailed description of the steps of performing flow calculation by thecomputer is given below.

Step S1401 of specifying a model gas system is given, for example:

Referring again to FIG. 6 showing a schematic view in section of themodel gas in the stable between parallel infinite plates in the openchannel at the uniform temperature, a case when both the lower plate 616and the upper plate 615 plates are at rest. The model gas is Newtonianwith constant viscosity. The pressure gradient is constant along withthe model gas flow in the x-direction. The main control volume hasdimensions Δx-length, Δy-thickness, and Δz-width. The pressure appliedto the main control volume at x and x+Δx are P_(x) and −P_(x+Δx),respectively. The velocity profile induced in the model gas due to themotion of the plates and pressure gradient is to be determined.

In a stable incompressible model gas system at the uniform temperatureand with a lack of internal processes, (1) each particle conservesproperties acquired by it at a point and time of the preceding collisionand deliver that properties unchanged in another location upon itssecond collision with another particle, (2) particle density n=constant,(2) temperature T and the magnitude of thermal velocity v_(T) areconstant. (3) We also accept, from Simulation 1, that the mass flowvelocity vector along y-direction u=0. For simplicity, the particles areconsidered to have a unit mass.

These parameters describing the model gas system are specified:

kind of the model gas including the mass of particles and an effectivecross-section of the collisions;

temperature T, pressure P;

the distance between parallel plates, H;

The pressure gradient along the channel,

$\frac{dP}{dx},$

and

allowable error of computation, Max Er.

Step S1402 of forming u_(x)-momentum balance equation includes, forexample, steps of:

(1) calculating the net rate of property influx, B_(in) ^(Ψ0_F), whichis formed by converging ballistic particles chosen from the group ofparticles originated from preceding collisions within space occupied bythe model gas (ballistic trajectories 601 and 603 in FIG. 6), fromEquation (124) given below, obtained by substitution of Equation (58)and assigning λ=0, Z_(V)=constant, u=0, B_(in) ^(Ψλ_F)≡B_(in) ^(Ψ0_F),Ψ_(in)=Ψ_(λ)≡Ψ(y′), and v(t_(i)′,y′,t,y)≡v_(T)n_(i) in Equation (38):

$\begin{matrix}{{{B_{in}^{\Psi \; 0\; \_ \; F}\left( {y,t} \right)} = {{- \frac{1}{2}}Z_{V}\frac{\partial}{\partial y}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}n_{i}{\Psi \left( y^{\prime} \right)}d\; y^{\prime}}}}},{where}} & (124) \\{{n_{i} = \frac{y - y^{\prime}}{{y - y^{\prime}}}},} & (125)\end{matrix}$

(2) representing a general integral form of the property balance byapplying Equation (126) given below, which is obtained by substitutionof Equations (124), (60), (63) in Equation (85):

$\begin{matrix}{{{{{- \frac{1}{2}}Z_{V}\frac{\partial}{\partial y}{\int_{0}^{H}{{Q_{i}\left( {y,y^{\prime}} \right)}n_{i}{\Psi \left( y^{\prime} \right)}{dy}^{\prime}}}} + {{\exp \left( {- {P_{c}\left( {H - y} \right)}} \right)}P_{c}\frac{1 - \sigma}{1 - {\left( {1 - \sigma} \right)^{2}{\exp \left( {{- 2}P_{c}H} \right)}}}\frac{1}{2}Z_{V}{\int_{0}^{H}{{\Psi \left( y^{\prime} \right)}{\exp \left( {{- P_{c}}y^{\prime}} \right)}{dy}^{\prime}}}} + {{\exp \left( {{- P_{c}}y} \right)}P_{c}\frac{1 - \sigma}{1 - {\left( {1 - \sigma} \right)^{2}{\exp \left( {{- 2}P_{c}H} \right)}}}\frac{1}{2}{\exp \left( {{- P_{c}}H} \right)}Z_{V}{\int_{0}^{H}{{\Psi \left( y^{\prime} \right)}{\exp \left( {P_{c}y^{\prime}} \right)}{dy}^{\prime}}}} + {{\overset{.}{Q}}_{\Psi}(y)}} = {Z_{V}{u_{x}(y)}}};} & (126)\end{matrix}$

(3) forming u_(x)-momentum balance equation by forming an integral formof the u_(x) momentum balance equation in a stable model gas flowconfined between two parallel plates with mixed diffuse and specularscatterings from the plates being at rest includes applying Equation(127) given below, which is obtained by substitution in Equation (126)of Ψ≡u_(x),

${{\overset{.}{Q}}_{\Psi} = {- \frac{dP}{dx}}},$

J₂₁ ¹ as of Equation (61), J₁₂ ² as of Equation (62),

as of Equation (54), and following differentiation and normalization byz_(V)≠0:

$\begin{matrix}{u_{x} = {{{- \frac{1}{z_{V}}}\frac{dP}{dx}} + {\frac{1}{2}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{u_{x}\left( y^{\prime} \right)}{dy}^{\prime}}}} + {\frac{1}{2}P_{c}\frac{\begin{matrix}{\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}} \\\left\lbrack {{\exp \left( {- {P_{c}\left( {H - y} \right)}} \right)} + {\exp \left( {{- P_{c}}y} \right)}} \right\rbrack\end{matrix}}{1 - {\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}}}{\int_{0}^{H}{{u_{x}\left( y^{\prime} \right)}{\exp \left( {P_{c}y^{\prime}} \right)}{{dy}^{\prime}.}}}}}} & (127)\end{matrix}$

Step S1403 of formulating solutions for the first approximationincludes, for example, steps of:

(1) assigning in the right-hand term of Equation (127)

u _(x)=0,  (128)

which is the velocity of each particle surrounding point y to be equalzero,

so Equation (127) is reduced to

$\begin{matrix}{{{u_{x\; 1}(y)} = {g_{0} = {{- \frac{1}{z_{V}}}\frac{dP}{dx}}}},} & (129)\end{matrix}$

which is a solution for the first order approximation;

(2) modification of the u_(x)-momentum balance equation, which is shownby Equation (127) in such form that u_(x) in the left-hand of theequation defines k^(th) approximation, u_(xk), and u_(x)(t_(i)′,y′) inthe right-hand of the equation defines (k−1)^(th) approximation,u_(x(k−1)), so Equation (127) is modified to

$\begin{matrix}{\mspace{79mu} {{u_{xk} = {{g_{({k - 1})}(y)} + g_{0}}},{where}}} & (130) \\{\mspace{85mu} {{g_{0} = {{- \frac{1}{z_{V}}}\frac{dP}{dx}}},{and}}} & (131) \\{{g_{{({k - 1})} +}(y)} = {{{+ \frac{1}{2}}P_{c}{\int_{0}^{H}{{\exp \left( {{- P_{c}}{{y - y^{\prime}}}} \right)}{u_{x{({k - 1})}}\left( y^{\prime} \right)}{dy}^{\prime}}}} + {\frac{1}{2}P_{c}\frac{\begin{matrix}{\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}} \\\left\lbrack {{\exp \left( {- {P_{c}\left( {H - y} \right)}} \right)} + {\exp \left( {{- P_{c}}y} \right)}} \right\rbrack\end{matrix}}{1 - {\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}}}{\int_{0}^{H}{{u_{x{({k - 1})}}\left( y^{\prime} \right)}{\exp \left( {P_{c}y^{\prime}} \right)}{{dy}^{\prime}.}}}}}} & (132)\end{matrix}$

Here g₀(y) is a well-defined function.

Step S1404 of formulating solutions for the k^(th) approximationincludes, for example, steps of:

(1) presenting the approximated solution for the n^(th) approximation,u_(xn), for each point of the space occupied by the model gas byEquations (133) given below:

u _(xk)(y)=u _(x1)(y)+u _(ak)(y),  (133)

where u_(a_n)(y) is the adjustment to velocity u_(x(n−1))(y) associatedthe velocities of particles surrounding point y, u_(x(n−1)+) (y′), whichis expressed as follows

$\begin{matrix}{{u_{ak}(y)} = {{\frac{1}{2}P_{c}{\int_{0}^{y}{{\exp \left( {- {P_{c}\left( {y - y^{\prime}} \right)}} \right)}{u_{{x{({k - 1})}} +}\left( y^{\prime} \right)}{dy}^{\prime}}}} + {\frac{1}{2}P_{c}{\int_{y}^{H}{{\exp \left( {P_{c}\left( {y - y^{\prime}} \right)} \right)}{u_{{x{({k - 1})}} +}\left( y^{\prime} \right)}{dy}^{\prime}}}} + {\frac{1}{2}P_{c}\frac{\begin{matrix}{\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}} \\\left\lbrack {{\exp \left( {- {P_{c}\left( {H - y} \right)}} \right)} + {\exp \left( {{- P_{c}}y} \right)}} \right\rbrack\end{matrix}}{1 - {\left( {1 - \sigma} \right){\exp \left( {{- P_{c}}H} \right)}}}{\int_{0}^{H}{{u_{{x{({k - 1})}} +}\left( y^{\prime} \right)}{\exp \left( {P_{c}y^{\prime}} \right)}{dy}^{\prime}}}}}} & (134)\end{matrix}$

The magnitude of u_(x)-momentum in a specific point in the model gasspace results from the impact of both the rate of momentum influx formedby particles being directly affected by the pressure force and the rateof momentum influx formed by particles that have acquired the propertyfrom each of their preceding collision being within the space occupiedby the model gas.

Step S1405 of establishing discretization parameters includes, forexample, steps of:

specifying discretization parameters including the dimension of thesystem for calculation [0,H] and the number of space intervals S; and

dividing the space intervals [0,H] into s equal parts of length l, sothat, for y_(i) ≡y′=sl>y, the range is divided as y=y_(s)<y_(s+1)< . . .<y_(s)=H and, for y_(s) ≡y′=sl<y, the range is divided as 0=y₀<y₁< . . .<y_(s)=y.

Step S1406 of providing an error estimation rule includes, for example,the step of providing the error estimation rule for k^(th) cycle of thesequential approximation from Equation (135) given below:

$\begin{matrix}{{Er} = {\sqrt{\frac{1}{s}{\sum_{1}^{k}\left( \frac{{u_{ak}\left( y_{s} \right)} - {u_{a{({k - 1})}}\left( y_{s} \right)}}{u_{xs}\left( y_{s} \right)} \right)^{2}}}.}} & (135)\end{matrix}$

Step S1407 of computing and storing the values u_(x1s)(y_(s)) forcorresponding y_(s) is done by applying Equation (129) given above;

Step 1408 of executing computation for each cycle k of the sequentialapproximations includes, for example, steps of:

assigning u_(xks+)=u_(x1s) (S1408A);

storing the values of u_(xks+) for corresponding Ys (S1408B);

computing u_(aks) ≡u_(ak)(y_(s)) for corresponding Ys according toEquation (134) (S1408C);

computing u_(xks)=u_(aks)+u_(x1s) according to Equation (133) andstoring the values of u_(xks) for corresponding y_(s) (S1408D);

computing Er according to Equation (135) given above (S1408E);

defining a condition for ending approximations (S1408F),

where If Er>Max_Er, assigning u_(xks+)=u_(xks), and go to (S1408B) and

if ≤Max_Er, the cycle of the sequential approximations is stopped, thencopying u_(xks) for corresponding y_(s), and

displaying the results of computation (S1409).

The sign + in the subscript indicates that the function is used for thefollowing integration with y′ as an argument.

FIG. 16 is a chart showing the gas flow velocity profiles across betweenthe plates at different numbers of the sequential approximations withmixed diffuse and specular scattering from the plates, which areobtained by computer simulations performed by the steps listed above,for stable conditions involving incompressible gas at the uniformtemperature. The computation involves the case of mixed diffuse andspecular scattering from the plates. Also, a comparative analyticalvelocity profile is obtained independently by using the analytical toolsdisclosed in the specification. As shown, the result of 120 sequentialapproximations computed according to Simulation 2 coincides with theanalytical velocity profile. The numerical calculations were made forthe nitrogen gas flow in a channel confined by parallel plates, whichare characterized by the property accommodation coefficient of 0.7 andspaced at the distance of 0.1 meters, at the uniform temperature of 300Kand pressure of 1 Pa. The simulations demonstrate significantimprovement from the prior art, which in addition to the capability tomodel and compute flows in rarefied gases also consists in the abilityto interpret the effect of the property/momentum accommodationcoefficient quantitatively.

In the above, Simulation 1 and Simulation 2 and comparative flowcalculations by a method of differentiation show an implementation ofembodiments of the which can be used for simulating transport processesin model gas systems.

The simulations set forth above are provided to give those of ordinaryskill in the art a complete disclosure and description of how to makeand use the embodiments of the methods and systems for predicting theevolution of the model gas fluid, which is based on the dynamicevolution of the model gas system in some period of time preceding theinitial time and do not limit the present invention. Modifications ofthe above-described modes for carrying out the disclosure may be used bypersons of skill in the art and are within the scope of these claims.

2.3 Example 3: A Special Purpose Computer for Implementing the Methodfor Simulating Fluid Flow

In some embodiments, the calculations are performed by using a computerprogram that directs a general purpose computer system to perform theactions of the above-described methods and analytical tools with, so thegeneral purpose computer is adapted by the computer program to operateas a special purpose computer for simulating the fluid flow. FIG. 17 isa schematic representation of the special purpose computer forimplementing the method for simulating the fluid flow. The specialpurpose computer 1700, for example, includes: a processor 1705, whichcan be any suitable computer processor; a non-transitory computerreadable memory coupled to the processor 1704, which can be any suitablecomputer readable and programmable memory; one or more user inputdevices such as a keyboard 1702 and/or a pointing device such as a mouse1703 operatively coupled to the processor via any suitable user I/Ointerface 1701; a display 1706, which can be any suitable computerdisplay, operatively coupled to the processor; the computer program,which is configured to use the fluid model and the model of propertybalance, may be stored in the non-transitory computer readable memory1704 and executable by the processor 1705. The computer program isoperable to run on the computer to execute a sequence of instructionsincluding:

setting parameters specifying geometry model of a fluid system, which isapproximated to a model gas system, including configurations,dimensions, materials, transferring temperature, pressure, and temporalcharacteristics involving a pre-initial time and an initial time;

transferring into memory a set of discretization parameters, wherein theset of discretization parameters comprises a set of points including aset of space points in space occupied by a model gas and a set ofsurface points on one or more of gas-solid interfaces being in contactwith the model gas, the set of points in which the model gas system isidentified by defining model gas material, pressure, temperature,external forces, and other pre-established properties comprising one ormore of mass, momentum, and energy during a time period from thepre-initial time until the initial time;

creating databases linking a selected non-moving point of the set ofspace points to each of the other points of the set of points;

computing, for each of links, a corresponding converging ballistictrajectory of a particle having a starting point in one point of the setof points and targeting the selected non-moving point at a given time;

computing for each of corresponding converging ballistic trajectories astarting time;

computing a local initial time in each of the set of points;

specifying property value obtained by each of converging ballisticparticles at the starting time, wherein, when the starting time is lessthan the local initial time, then obtained property value is defined bya pre-established property content;

computing probability of traveling along each of correspondingconverging ballistic trajectories;

computing a net rate of property influx per unit volume, which is formedby way of a plurality of converging ballistic particles, in each of aset of non-moving points at the given time from the model gas system;

computing for each link a corresponding diverging ballistic trajectoryof a diverging ballistic particle escaping the selected non-moving pointat the given time; identifying property obtained by each of divergingballistic particles at the given time;

computing probability of traveling along each of corresponding divergingballistic trajectories;

computing a net rate of property efflux per unit volume, which is formedby way of a plurality of diverging ballistic particles, from each of theset of non-moving points at the given time into the model gas;

establishing a property balance equation by formulating balance betweenthe net rate of property influx and the net rate of property efflux ineach of the set of non-moving points at the given time, therebyformulating a system of mutually dependent integral property balanceequations; and

resolving the system of mutually dependent integral property balanceequations and obtaining property value comprising flow velocity value ineach of the set of space points, thereby reliably predicting gas flow,including rarefied gas flows.

The method, for example, is implemented in software in any suitablesoftware language, including C++, Visual C, Visual C++, Visual Basic,Java, C, and FORTRAN. The software program may be stored on any computerreadable medium, including a floppy disk, a hard disk drive, amagneto-optical disk drive, CD-ROM, magnetic tape or any other ofseveral non-volatile storage devices well known to those skilled in theart for storing a computer software program in such a form that theinstructions may be retrieved and executed by a processor.

INDUSTRIAL APPLICABILITY

The most significant improvements in CFD technology by the methoddisclosed above will advance the following applications and industries:

1. Efficient and accurate modeling of transport processes in rarefiedgases according to the inventive approach will benefit airspaceindustries, astrophysics, microelectronics and nanoelectronicsprocessing, and plasma processing.

2. Enhancement in time resolution of computations for unstable flowsaccording to the inventive approach will advance applications forturbulence, shockwave, and detonation analysis, which in turn will boostsupersonic aircraft and spacecraft development.

3. Innovative quantitative interpretation of the aging factor accordingto the inventive method will enhance applications for radiation andnuclear physics, physics of chemical reactions, and plasma physicsprocessing.

4. Finally, an improvement associated with the disclosed method oftaking into account a pre-established or known dynamic history of thefluid system during a pre-initial period will benefit any application byproviding increased certainty and accuracy of prediction of the fluidsystem evolution.

A new generation of Computational Fluid Dynamics software products,which are based on the disclosed analytical tools and method, areanticipated having capability to modeling gases from the continuum flowregime (Kn<0.01) to the rarefied and free molecular flow regime (Kn>10),considerably higher accuracy of prediction, and lower computation cost.

It is to be understood that the present invention is not limited to theembodiments described above, but encompasses any and all embodimentswithin the scope of the following claims.

1. A computer implemented method for modeling transport processes influids comprising: approximating fluid flow in a fluid system as a flowof a model gas in a model gas system being identical to the fluidsystem; assuming that each of a plurality of particles, which composethe model gas, travels with a probability between any points in a spaceoccupied by the model gas by following a ballistic trajectory governedby a law of motion in free space, the ballistic trajectory having astarting point in one of a plurality of points of original collisionsand an ending point in one of a plurality of points of endingcollisions; treating each of the plurality of points of originalcollisions as a point source; treating each of the plurality of pointsof ending collisions as a point sink; treating each of a plurality ofballistic particles moving from the point source to the point sink as aproperty carrier created in the point source at a time of an originalcollision by obtaining one or more of properties comprising one or moreof mass, momentum, and energy of specific values being intrinsic to thepoint source, and ended in the point sink at the time of an endingcollision by transferring one or more of properties of specific valuesin the point sink; forming an integral property balance equation foreach of one or more properties being transported by the plurality ofballistic particles in each of a plurality of non-moving points at agiven time, wherein the plurality of ballistic particles in each of theplurality of non-moving points at the given time comprises a pluralityof converging ballistic particles and a plurality of diverging ballisticparticles, and wherein each of the plurality of non-moving points istreated as the point sink for each of the plurality of convergingballistic particles and as the point source for each of the plurality ofdiverging ballistic particles; and
 1. A computer implemented method formodeling transport processes in fluids comprising: using a fluid model,wherein the fluid model comprises treating a fluid flow including a flowof rarefied gases in a fluid system as a flow of a model gas in a modelgas system being identical to the fluid system, wherein the model gas iscomposed of a plurality of particles including molecules, which moverandomly and interact by collisions, wherein each of the plurality ofparticles is assigned to travel with a probability between any points ina space occupied by the model gas by following a ballistic trajectorygoverned by a law of motion, the ballistic trajectory having a startingpoint in a point of original collision and an ending point in a point ofending collision, wherein each of the plurality of points of originalcollision is treated as a point source, wherein each of the plurality ofpoints of ending collision is treated as a point sink, and wherein eachof a plurality of ballistic particles moving from the point source tothe point sink is treated as a property carrier created in the pointsource at a time of an original collision by obtaining one or more ofproperties comprising one or more of mass, momentum, and energy ofspecific values being intrinsic to the point source, and ended in thepoint sink at a time of an ending collision by transferring one or moreof properties of specific values in the point sink; and using a computerfor simulating the fluid flow, wherein the computer for simulating thefluid flow comprises a computer readable medium for storing dataassociated with a computer program, and a processor being connected foraccessing the computer program and operable to run the computer, whereinthe computer program is configured to use the fluid model and a model ofproperty balance, wherein the model of property balance in a combinationthe fluid model transforms parameters characterizing motion of aplurality of ballistic particles into parameters characterizing thefluid flow, wherein the model of property balance comprises establishingan integral property balance equation for a property being transportedby a plurality of converging ballistic particles and a plurality ofdiverging ballistic particles in each of a plurality of non-movingpoints in the space occupied by the model gas at a given time, andwherein the computer program is operable to run on the computer forsimulating the fluid flow to compute the flow of the model gas in eachof the plurality of non-moving points at the given time and to displayresults of computing, thereby reliably predicting the fluid flowincluding the flow of rarefied gases.
 2. The method of claim 1, whereina value of property delivered in the point sink by each of the pluralityof converging ballistic particles is evaluated regarding whether thevalue of property is conserved, whether the value of property is changedbecause of aging, and whether the value of property is modified becauseof interaction with an external field during a ballistic traveling time.3. The method of claim 1, wherein each of a plurality of collisions on agas-solid interface of the model gas system, which has resulted indiffuse particle scattering from the gas-solid interface, is treated asan act of interaction involving a property transfer from the gas-solidinterface to a scattered particle, wherein each of a plurality of pointsof diffuse particle scattering on the gas-solid interface is treated asa heterogeneous point source for each of a plurality of scatteredparticles, wherein velocity of each of a plurality of heterogeneouspoint sources on the gas-solid interface is assigned to be equal to thevelocity of the gas-solid interface in each of a plurality ofcorresponding points of diffuse particle scattering at time of diffuseparticle scattering, and wherein point source strength of each of theplurality of heterogeneous point sources on the gas-solid interface isassigned to be directly proportional to a property accommodationcoefficient in each of the plurality of corresponding points of diffuseparticle scattering at the time of diffuse particle scattering, theproperty accommodation coefficient which is the probability, for anincident particle, to accommodate one or more of properties intrinsic tothe gas-solid interface and to scatter back in the model gas, theproperty accommodation coefficient being in a range from zero to one. 4.The method of claim 1, further comprising steps of: establishing aninitial time, wherein the initial time is treated as the time of aspecific modification of the model gas system in a specific locationcomprising specific modification of one or more of temperature, velocityof a gas-solid interface, and application of an external field that maymodify a property value; specifying the given time, wherein the giventime is greater than or equal to the initial time; specifying the modelgas system, wherein the model gas system is specified by defining modelgas properties comprising a gas material, pressure, temperature and byestablishing a geometry model, wherein establishing the geometry modelcomprises setting geometry and boundary conditions during a period froma pre-initial time until the given time; establishing discretizationparameters comprising a set of discretization points in the spaceoccupied by the model gas, wherein each of the set of discretizationpoints is assigned to a corresponding point of the plurality ofnon-moving points; calculating, for each of the set of discretizationpoints in the space occupied by the model gas, a local initial time,wherein the local initial time is greater than or equal to the initialtime; obtaining, for each of the set of discretization points in thespace occupied by the model gas at a given advanced time, which is aheadof the local initial time, a past property value; formulating, for eachof the set of discretization points at the given time, an approximatedintegral form of property balance, wherein the approximated integralform of property balance is established in a way that, in a selectedpoint of the set of discretization points at the given time, a net rateof property influx per unit volume from the model gas system, which isformed by a set of converging particles and calculated by summing a rateof property influx from the set of discretization points surrounding theselected point is equated to a net rate of property efflux per unitvolume from the selected point of the set of discretization points,separating the plurality of converging ballistic particles on a firstplurality of converging ballistic particles and on a second plurality ofconverging ballistic particles, where each of the first plurality ofconverging ballistic particles delivers one or more of past propertyvalues from the model gas system, where each of the second plurality ofconverging ballistic particles delivers one or more of present propertyvalues from the model gas system, where each of past property values isobtained from one of the plurality of points of original collisions atthe time ahead of the local initial time, where each of past propertyvalues is well defined in each of the plurality of points of originalcollision at time of original collision, where each of present propertyvalues is obtained from one of the plurality of points of originalcollisions at the time which follows the local initial time, and where apresent property value in each of the plurality of points of theoriginal collision is treated as an unknown property value; andcomputing the flow of the model gas by considering an impact of aplurality of initial converging ballistic particles and taking intoaccount an effect of mutual dependence of property exchange betweenparticles interacting by collisions, which comprises a sequentialapproximation method.
 5. The method of claim 1, wherein the integralproperty balance equation is established in a way that, in a generalpoint of the plurality of non-moving points in the space occupied by themodel gas at the given time, a net rate of property influx per unitvolume, which is formed by the plurality of converging ballisticparticles, each of the plurality of converging ballistic particles isselected from the plurality of ballistic particles by the ballistictrajectory having the ending point in the general point at the giventime, is equated to sum a temporal rate of a property change per unitvolume and a net rate of property efflux per unit volume, which isformed by the plurality of diverging ballistic particles, each of theplurality of diverging ballistic particles is selected from theplurality of ballistic particles by the ballistic trajectory having thestarting point in the general point at the given time.
 6. A computerimplemented method for modeling transport processes in fluidscomprising: using a fluid model, wherein using the fluid model comprisestreating a fluid flow including a flow of rarefied gases in a fluidsystem as a flow of a model gas in a model gas system being identical tothe fluid system, wherein the model gas is composed of a plurality ofparticles, including molecules, which move randomly and interact bycollisions, wherein each of the plurality of particles is assigned totravel with a probability between any of two points in a space occupiedby the model gas by following a ballistic trajectory governed by a lawof motion, the ballistic trajectory having a starting point in a pointof original collision and an ending point in a point of endingcollision, and wherein each of a plurality of ballistic particlestransports a combination of one or more of properties comprising one ormore of mass, momentum, and energy from the starting point to the endingpoint; using a model of property balance, wherein the model of propertybalance in a combination the fluid model transforms parameterscharacterizing motion of a plurality of ballistic particles intoparameters characterizing the fluid flow, wherein the model of propertybalance comprises: specifying a geometry model, defining a net rate ofproperty influx per unit volume in a general non-moving point in thespace occupied by the model gas at a given time, the net rate ofproperty influx per unit volume, which is formed by a plurality ofconverging ballistic particles, each of the plurality of convergingballistic particles is selected from the plurality of ballisticparticles by the ballistic trajectory having the ending point in thegeneral non-moving point at the given time, defining a net rate ofproperty efflux per unit volume from the general non-moving point at thegiven time, the net rate of property efflux per unit volume, which isformed by a plurality of diverging ballistic particles, each of theplurality of diverging ballistic particles is selected from theplurality of ballistic particles by the ballistic trajectory having thestarting point in the general non-moving point at the given time, andestablishing an integral property balance equation for each of one ormore properties being transported by the plurality of ballisticparticles in the general non-moving point; and using a computer forsimulating the fluid flow, wherein the computer for simulating the fluidflow comprises a computer readable medium for storing data associatedwith a computer program, and a processor for executing a sequence ofinstructions from the computer program, wherein the computer program isconfigured to use the fluid model and the model of property balance, andwherein the computer program is operable to run on the computer forsimulating the fluid flow to compute the flow of the model gas in everygeneral point of a plurality of non-moving points in the space occupiedby the model gas at the given time and to display results of computing,thereby reliably predicting the fluid flow including the flow ofrarefied gases.
 7. The method of claim 6, wherein the step of specifyingthe geometry model further comprises: establishing geometry and boundaryconditions, wherein geometry and boundary conditions are establishedduring a period from a pre-initial time until the given time; andestablishing a dynamic history of the model gas, wherein the dynamichistory of the model gas, which is in form of one or more of instantdistribution patterns of model gas properties within the model gassystem at a set of different sequential moments in time, is establishedduring a pre-initial period, and wherein the pre-initial period is setto be not less than a reliable period before an initial time, whereinthe initial time is treated as a time before which the fluid system isdescribed,
 8. The method of claim 7, wherein the reliable period beforethe initial time is estimated as${\phi_{i\; 0}^{rel} \cong {- \frac{\ln ({MAV})}{\left\lbrack {P_{c}v_{rel}} \right\rbrack_{m}}}},$where φ_(i0) ^(rel) is the reliable period before the initial time, MAVis a minimum absolute value corresponding to a precision formatsupported by a computer hardware, P_(c) is a number of particles withina collision tube of a unit length placed in the model gas, v_(rel) is anaverage magnitude of velocity of a traveling particle with respect to anearby particle, and [P_(c)v_(rel)]_(max) is a highest value of productof P_(c)v_(rel) in the model gas system.
 9. The method of claim 6,wherein the plurality of converging ballistic particles comprises afirst plurality of converging ballistic particles, where the ballistictrajectory of one of the first plurality of converging ballisticparticles has the starting point in the space occupied by the model gas,wherein the step of defining the net rate of property influx per unitvolume from the model gas system further comprises a step of defining anet rate of total property influx per unit volume from a surroundingmodel gas, wherein the net rate of total property influx per unit volumeis formed by a first plurality of converging ballistic particles, eachof the first plurality of converging ballistic particles is selectedfrom the plurality of converging ballistic particles by the ballistictrajectory having the starting point in one of a plurality of points oforiginal collisions within the space occupied by the model gas, andwherein each of the plurality of points of original collisions istreated as a point source.
 10. The method of claim 9, wherein the stepof defining the net rate of total property influx per unit volume fromthe surrounding model gas comprises: identifying, for each of the firstplurality of converging ballistic particles, the ballistic trajectory;defining the probability of traveling along the ballistic trajectory;defining a net rate of particle efflux per unit volume from one of theplurality of points of an original collision at time of the originalcollision, wherein the time of the original collision for each of theplurality of converging ballistic particles is identified by theballistic trajectory, and wherein velocity of the point source for eachof the first plurality of converging ballistic particles is assigned tobe equal to mass flow velocity of the model gas in a corresponding pointof the original collision at the time of the original collision;defining a net rate of property flux per unit area in the general pointat the given time from one of a plurality of point sources; defining anet rate of total property flux per unit area in the general point atthe given time from the plurality of point sources; and defining the netrate of total property influx per unit volume in the general point atthe given time by applying divergence operator to the net rate of totalproperty flux.
 11. The method of claim 10, wherein the step of theidentifying the ballistic trajectory comprises: determining one or moreof characteristics of the ballistic trajectory comprising the time ofthe original collision, an appropriate instant unit vector defining theballistic trajectory in the starting point, and a velocity vector in theending point of the ballistic trajectory at the given time; and defininga size of an expansion zone.
 12. The method of claim 10, wherein thestep of defining the net rate of total property influx per unit volumein a one-dimensional configuration further comprises a step offormulating the net rate of total property influx as${{B_{in}^{\Psi \; \lambda \; \_ \; F}\left( {y,t} \right)} = {{- \frac{1}{2}}\frac{\partial}{\partial y}{\int_{{yb}\; 2}^{{yb}\; 1}{{Q_{i}\left( {t,t_{i}^{\prime}} \right)}{Z_{V}\left( {t_{i}^{\prime},y^{\prime}} \right)}\frac{v\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{v_{T}\left( {t_{i}^{\prime},y^{\prime}} \right)}{\Psi_{in}\left( {t_{i}^{\prime},y^{\prime},t,y} \right)}{dy}^{\prime}}}}},$wherein $\frac{\partial}{\partial y}$ is the divergence operator inone-dimensional configuration, B_(in) ^(Ψλ_F) is the net rate ofproperty influx per unit volume, t is the given time, t_(i)′ is the timeof the original collision, y is position of the ending point in one ofthe plurality of points of ending collisions, y′ is position of thestarting point in one of the plurality of points of original collisions,v_(T)(t_(i)′,y′) is average magnitude of thermal velocity in thestarting point at time of original collision, Z_(V)(t_(i)′,y′) is rateof collisions per unit volume in the starting point at the time oforiginal collision, v(t_(i)′,y′,t,y) is a velocity vector in the endingpoint at the given time, Q_(i)(t,t_(i)′) is the probability of travelingalong the ballistic trajectory, Ψ_(in)(t_(i)′,y′,t,y) is a propertyvalue delivered in the ending point at the given time by one of thefirst plurality of converging ballistic particles, yb1 is a lower limitof integration over volume of space occupied by the model gas, yb2 is anupper limit of integration over volume of space occupied by the modelgas, and Z_(V)(t_(i),y′) is calculated asZ _(v)(t _(i) ′,y′)=½P _(c)(t _(i) ′,y′)v _(rel)(t _(i) ′,y′)n(t _(i)′,y′), where n(t_(i)′,y′) is particle density, P_(c)(t_(i)′,y′) is anaverage number of collisions per unit length, and v_(rel)(t_(i)′,y′) isan average magnitude of relative velocity.
 13. The method of claim 12,wherein the probability of traveling along the ballistic trajectory inone-dimensional configuration is calculated asQ _(i)(t,t _(i)′)=exp(−∫_(t) _(i) _(′) ^(t) P _(c)({tilde over(y)}({tilde over (t)}))v _(rel)({tilde over (y)}({tilde over(t)}))d{tilde over (t)}), where {tilde over (t)} is a parametric time ina range of t_(i)′<{tilde over (t)}≤t, {tilde over (y)}({tilde over (t)})is a trajectory point of the ballistic trajectory at parametric time{tilde over (t)}, v_(rel)({tilde over (y)}({tilde over (t)})) is theaverage magnitude of relative velocity in the trajectory point, andP_(c)({tilde over (y)}({tilde over (t)})) is average number ofcollisions per unit length in the trajectory point.
 14. The method ofclaim 13, wherein the average magnitude of relative velocity in thetrajectory point is calculated by steps of: defining an instantmagnitude of velocity of one of the plurality of converging ballisticparticles in the trajectory point with respect to a nearby particle inthe trajectory point; and averaging the instant magnitude of velocityover a plurality of directions of a thermal velocity component of anearby particle.
 15. The method of claim 6, wherein a property valuedelivered in the ending point by each of the plurality of convergingballistic particles is calculated asΨ_(in)=Ψ_(λ)+Ψ_(g), where Ψ_(in) is the property value delivered in theending point, Ψ_(g) is a field function characterizing property value ofknown measure, which is changed during traveling time, φ, because ofinteraction with external field comprising modification of momentum ingravitation field, Ψ_(λ) is an age function characterizing the propertyvalue in the ending point of the ballistic trajectory, which iscalculated asΨ_(λ)=Ψ₀ e ^(−λφ), where Δ is an aging coefficient, which is in rangefrom negative to positive values including zero value, and where Ψ₀ is astarting function characterizing the property value carried by each ofthe plurality of ballistic particles in the starting point of theballistic trajectory at a starting time, t_(i), which is calculated asΨ₀(t _(i) ′,y′)=Ψ⁰⁻[t _(i) <t _(i0)]+Ψ[t _(i) ≥t _(i0)], where t_(i0) isa local initial time, each of expressions with square brackets isIverson bracket, which converts a statement in brackets into a number,the number is one if the statement is satisfied, and the number is zerootherwise, Ψ⁰⁻ is a pre-established property value in each of theplurality of points in space of the model gas system at the startingtime ahead of the local initial time, and Ψ is a present property valuein each of the plurality of points in space of the model gas system atthe starting time greater than or equal to the local initial time. 16.(canceled)
 17. The method of claim 6, wherein the plurality ofconverging ballistic particles comprises a second plurality ofconverging ballistic particles, where the ballistic trajectory of one ofthe second plurality of converging ballistic particles has the startingpoint on a gas-solid interface, wherein the step of defining the netrate of property influx per unit volume from the model gas systemfurther comprises a step of defining a net rate of total property influxper unit volume from a gas-solid interface exhibiting mixed diffuse andspecular particle scattering, wherein the net rate of total propertyinflux per unit volume from the gas-solid interface is formed by way ofa second plurality of converging ballistic particles, each of the secondplurality of converging ballistic particles is selected from theplurality of converging ballistic particles by the ballistic trajectoryhaving the starting point in one of a plurality of points of originalcollisions on the gas-solid interface, and wherein each of the pluralityof points of original collisions resulted in diffuse particle scatteringfrom the gas-solid interface is treated as a heterogeneous point source.18. The method of claim 17, wherein the step of defining the net rate oftotal property influx per unit volume from the gas-solid interfacecomprises: identifying the ballistic trajectory, for each of the secondplurality of converging ballistic particles, wherein velocity of theheterogeneous point source for each of the second plurality ofconverging ballistic particles is assigned to be equal to the velocityof the gas-solid interface in a corresponding point of originalcollisions on the gas-solid interface at a time of diffuse particlescattering; defining the probability of traveling along the ballistictrajectory; defining a rate of collisions per unit area in one of aplurality of heterogeneous point sources on the gas-solid interface atthe time of diffuse particle scattering, wherein the time of diffuseparticle scattering for each of the second plurality of convergingballistic particles is identified by the ballistic trajectory; defininga net rate of property flux per unit area in the general point at thegiven time from one of the plurality of heterogeneous point sources;defining a net rate of total property flux in the general point at thegiven time from the plurality of heterogeneous point sources on one ormore of gas-solid interfaces; and defining the net rate of totalproperty influx per unit volume in the general point at the given timeby applying divergence operator to the net rate of total property flux.19. The method of claim 18, wherein the step of defining the rate ofcollisions comprises: defining a bulk component of an impingement rate,which is formed by particles chosen from the plurality of convergingballistic particles originated from the original collisions in the spaceoccupied by the model gas; defining a bulk-specular component of theimpingement rate, which is formed by particles chosen from the pluralityof converging ballistic particles originated from the originalcollisions in the space occupied by the model gas and having at leastthe last preceding specular scattering; defining a diffuse component ofthe impingement rate, which is formed by particles chosen from theplurality of converging ballistic particles originated from one or moreof heterogeneous point sources; defining a diffuse-specular component ofthe impingement rate, which is formed by particles chosen from theplurality of converging ballistic particles originated from originaldiffuse scatterings from one or more of heterogeneous point sources andhaving at least the last preceding specular scattering; and summing thebulk component, the bulk-specular component, the diffuse component, andthe diffuse-specular component of the impingement rate, wherein each ofa plurality of impinging particles is selected by an impingingtrajectory allowing to target the gas-solid interface in a specificlocation of the gas-solid interface at a specific time.
 20. The methodof claim 18, wherein the step of defining the net rate of total propertyinflux per unit volume further comprises a step of formulating the netrate of total property influx as${B_{in}^{\Psi \; \lambda \; \_ \; b} = {- {\frac{\partial}{\partial y}\left\lbrack {{Q_{ib}\left( {t,t_{ib}^{\prime}} \right)}\sigma \; {Z_{b}\left( {t_{ib}^{\prime},y_{b}} \right)}\frac{v_{b}\left( {t_{ib}^{\prime},y_{b},t,y} \right)}{v_{{Tb}\; 1}\left( {t_{ib}^{\prime},y_{b}} \right)}{\Psi_{inb}\left( {t_{ib}^{\prime},y_{b},\lambda,\phi_{{ib}\; 1}} \right)}} \right\rbrack}}},$wherein B_(in) ^(Ψλ_b) is the net rate of total property influx from thegas-solid interface in one-dimensional configuration, t is the giventime, t_(ib)′ is a time of scattering from the gas-solid interface, y isposition of the ending point, y_(b) is position of the starting point inone of the plurality of points of original collisions on the gas-solidinterface, v_(Tb)(t_(ib)′,y_(b)) is a magnitude of thermal velocityobtained in the starting point at the time of scattering from thegas-solid interface, Z_(b)(t_(ib)′,y_(b)) is the rate of collisions perunit area in the starting point at the time of scattering from thegas-solid interface, v_(b)(t_(ib)′,y_(b),t,y) is a velocity vector inthe ending point at the given time, Q_(ib)(t,t_(ib)) is the probabilityof traveling along the ballistic trajectory of the second plurality ofconverging ballistic trajectories, Ψ_(inb)(t_(ib)′,y′,λ,t,y) is aproperty value delivered in the ending point at the given time by one ofthe second plurality of converging ballistic particles, λ is agingcoefficient, and σ is a property accommodation coefficient.
 21. Themethod of claim 6, wherein the step of defining the net rate of propertyefflux per unit volume from the general point of the plurality ofnon-moving points at the given time in the model gas system comprises:identifying, for each of the plurality of diverging ballistic particles,the ballistic trajectory starting from the general point of theplurality of non-moving points at the given time and ending in one ofthe plurality of points of ending collisions surrounding the generalpoint; defining the probability of traveling along the ballistictrajectory; defining a vector field of a property flux per unit area inthe ending point of the ballistic trajectory of a diverging ballisticparticle; and defining the net rate of property efflux per unit volumeby applying a divergence operator to the vector field of the propertyflux per unit area in the ending point; and shrinking volume of spacesurrounding the general point to the general point.
 22. The method ofclaim 21, wherein the step of defining the net rate of property effluxper unit volume in one-dimensional configuration comprises a step offormulating the net rate of property efflux as${{B_{out}^{{\Psi\_}{FS}}\left( {t,y} \right)} = {\frac{1}{2}\left\{ {\frac{\partial}{\partial y}\left\lbrack {{{n\left( {t,y} \right)}Q} + {\left( {t_{a}^{\prime},t} \right)v} + {\left( {t,y,t_{a}^{\prime},y^{\prime}} \right){\Psi_{\lambda}\left( {t,y,t_{a}^{\prime},y^{\prime}} \right)}}} \right\rbrack} \right\}_{y^{\prime}\rightarrow y}}},$where $\frac{\partial}{\partial y}$ is the divergence operator inone-dimensional configuration, B_(out) ^(Ψ_FS) is the net rate ofproperty efflux per unit volume, t is the given time, y is position ofthe starting point of a plurality of ballistic trajectories of divergingballistic particles, y′ is position of the ending point of one of theplurality of ballistic trajectories of diverging ballistic particles,t_(a)′ is a time of positioning in the ending point, v₊(t,y,t_(a)′,y′)is a velocity vector of each of the plurality of diverging ballisticparticle at the time of positioning the ending point, Q₊(t_(a)′,t) isthe probability of traveling along the ballistic trajectory of one ofthe plurality of diverging ballistic particles, Ψ_(λ)(t,y,t_(a)′,y′) isproperty value carried by one of the plurality of diverging ballisticparticles at the time of positioning the ending point, and n(t,y) isparticle density in the starting point at the given time.
 23. A computerimplemented method for modeling transport processes in fluids being in astable condition comprising: using a fluid model, wherein using thefluid model comprises treating a fluid flow including a flow of rarefiedgases in a fluid system as a flow of a model gas in a model gas systembeing identical to the fluid system, wherein the model gas is composedof a plurality of particles, including molecules, which move randomlyand interact by collisions, wherein each of the plurality of particlesis assigned to travel with a probability between any of two points in aspace occupied by the model gas by following a ballistic trajectorygoverned by a law of motion, the ballistic trajectory having a startingpoint in a point of original collision and an ending point in a point ofending collision, and wherein each of a plurality of ballistic particlestransports a combination of one or more of properties comprising one ormore of mass, momentum, and energy from the starting point to the endingpoint; using a model of property balance, wherein the model of propertybalance in a combination the fluid model transforms parameterscharacterizing motion of a plurality of ballistic particles intoparameters characterizing the fluid flow, wherein the model of propertybalance comprises: specifying geometry and boundary conditions, defininga net rate of property influx per unit volume in a general non-movingpoint in the space occupied by the model gas, the net rate of propertyinflux per unit volume, which is formed by a plurality of convergingballistic particles, each of the plurality of converging ballisticparticles is selected from the plurality of ballistic particles by theballistic trajectory having the ending point in the general non-movingpoint, defining a net rate of property efflux per unit volume from thegeneral non-moving point, the net rate of property efflux per unitvolume, which is formed by a plurality of diverging ballistic particles,each of the plurality of diverging ballistic particles is selected fromthe plurality of ballistic particles by the ballistic trajectory havingthe starting point in the general non-moving point, and establishing anintegral property balance equation for each of one or more propertiesbeing transported by the plurality of ballistic particles in the generalnon-moving point; and using a computer for simulating the fluid flow,wherein the computer for simulating the fluid flow comprises a computerreadable medium for storing data associated with a computer program, anda processor for executing a sequence of instructions from the computerprogram, wherein the computer program is configured to use the fluidmodel and the model of property balance, and wherein the computerprogram is operable to run on the computer for simulating the fluid flowto compute the flow of the model gas in every general point of aplurality of non-moving points in the space occupied by the model gasand to display results of computing, thereby reliably predicting thefluid flow including the flow of rarefied gases.
 24. The method of claim23, further comprising: treating each point of a plurality of points inthe space occupied by the model gas as a point of collisions for each ofthe plurality of ballistic particles having the ballistic trajectorywith a same ending point; treating each of a plurality of points ofcollisions as either a point source for each of the plurality ofdiverging ballistic particles or a point sink for each of the pluralityof converging ballistic particles; and treating each of the plurality ofballistic particles moving from the point source to the point sink as aproperty carrier created in the point source by obtaining one or more ofproperties of specific values being intrinsic to the model gassurrounding the point source, and ended in the point sink bytransferring one or more of properties of specific values in the pointsink, wherein a value of property delivered by each of the plurality ofballistic particles converging in the point sink is evaluated regardingwhether the value of property is conserved, whether the value ofproperty is changed because of aging, and whether the value of propertyis modified because of interaction with an external field during aballistic traveling time.
 25. The method of claim 23, furthercomprising: treating each of a plurality of collisions on a gas-solidinterface of the model gas system, which has resulted in diffuseparticle scattering from the gas-solid interface, as an act ofinteraction involving a property transfer from the gas-solid interfaceto a scattered particle; and treating each of a plurality of points ofdiffuse particle scattering on the gas-solid interface as aheterogeneous point source for each of a plurality of scatteredparticles, wherein the gas-solid interface reveals mixed diffuse andspecular particle scatterings, wherein velocity of each of a pluralityof heterogeneous point sources on the gas-solid interface is assigned tobe equal to the velocity of the gas-solid interface in each of aplurality of corresponding points of diffuse particle scattering, andwherein point source strength of each of the plurality of heterogeneouspoint sources on the gas-solid interface is assigned to be directlyproportional to a property accommodation coefficient in each of theplurality of corresponding points of diffuse particle scattering, theproperty accommodation coefficient which is the probability, for anincident particle, to accommodate one or more of properties intrinsic tothe gas-solid interface and to scatter back in the model gas as adiffuse particle, the property accommodation coefficient being in arange from zero to one.
 26. (canceled)
 27. (canceled)
 28. A specialpurpose computer, comprising: a processor; a non-transitory computerreadable memory coupled to the processor; a user interface coupled tothe processor; a display coupled to the processor; and a program storedin the non-transitory computer readable memory and executable by theprocessor, wherein the program is operable to use a fluid model and amodel of property balance, wherein the fluid model comprises treating afluid flow including a flow of a rarefied gas as a flow of a model gasin a model gas system, the model gas in which particles includingmolecules are treated as randomly interacting by collisions ballisticparticles traveling with a probability between points in a spaceoccupied by the model gas by following a ballistic trajectory governedby a law of motion, wherein the model of property balance in acombination the fluid model transforms parameters characterizing motionof a plurality of ballistic particles into parameters characterizing thefluid flow, and wherein the program is operable when executed to:compute a trajectory and a velocity of a ballistic particle convergingin a general non-moving point at a given time; compute a probability ofballistic traveling along a converging ballistic trajectory; compute anet rate of property influx per unit volume in a general non-movingpoint in a space occupied by the model gas at the given time, which isformed by a flow of converging ballistic particles, each having apreceding collision in the space occupied by the model gas; compute thenet rate of property influx per unit volume in the general non-movingpoint of at the given time, which is formed by flow of convergingballistic particles, each having a preceding scattering from a gas-solidinterface; compute the trajectory and the velocity of the ballisticparticle diverging from the general non-moving point at the given time;compute the probability of ballistic traveling along a divergingballistic trajectory; and compute a net rate of property efflux per unitvolume, which is formed by diverging ballistic particles from thegeneral non-moving point at the given time, thereby upon formulating asystem of mutually dependent integral property balance equations in eachof a plurality of non-moving points at the given time and resolving thesystem of mutually dependent integral property balance equations,obtaining property value in each of the plurality of non-moving pointsat the given time, and thereby reliably predicting the fluid flowincluding the flow of the rarefied gas.